


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1996-06 


Frequency domain structural identification 


Johnson, Richard. 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/8428 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
. (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist L Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


nN i KNOX appointed — and published -- scholarly author. 
| inp 
, LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 





NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CALIFORNIA 





THESIS 


FREQUENCY DOMAIN STRUCTURAL 
IDENTIFICATION 


by 
Richard Johnson 


June 1996 


Thesis Advisor: Joshua H. Gordis 


Approved for public release; distribution is unlimited 





DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 





t 


> reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instruction, searching existing data sources, gathering and 
aining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, 
ling suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA} 
4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188) Washington DC 20503. 








AGENCY USE ONLY (Leave blank) , REPORT DATE 3: REPORT TYPE AND DATES COVERED 
June 1996 Master’s Thesis 
TITLE AND SUBTITLE 5: FUNDING NUMBERS 
EQUENCY DOMAIN STRUCTURAL IDENTIFICATION 
AUTHOR(S) Richard Johnson 





PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 
ral Postgraduate School REPORT NUMBER 
nterey CA 93943-5000 


SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the aad 





cy or position of the Department of Defense or the U.S. Government. . 
DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 
roved for public release; distribution is unlimited. 
ABSTRACT (maximum 200 words) 





; 
The Structural Synthesis Transformation is used to conduct structural system identification in the frequency domain. For spatially complete} 
s where each of the frequency response functions at every degree of freedom of each of the coordinates of the modeled system are available it} 
10wn that the theory exactly identifies all modeling errors. For spatially incomplete cases where the frequency response functions are! 
lable only at a proper subset of the degrees of freedom of the finite element model, single mode solutions are computed over intervals about] 
nodes of the experimental system using matrices and complex valued line integrals. Methods of forming multiple mode solutions from the] 
le mode solutions are explored. | 


/ 








SUBJECT TERMS Finite Element Method 15. NUMBER OF PAGES 
161 


16. PRICE CODE 













SECURITY CLASSIFICATION OF} 18. 
ORT 


lassified 


SECURITY  CLASSIFI- 
CATION OF THIS PAGE 


Unclassified 


19. 
CATION OF ABSTRACT 


Unclassified 












NSN 7540-01-28-5500 Standard Form 298 (Rev. 289) 


Prescribed by ANSI Std. 239-18 298-102 








Approved for public release; distribution is unlimited. 


FREQUENCY DOMAIN STRUCTURAL IDENTIFICATION 


Richard Johnson 
Lieutenant Commander, United States Navy 
B.S., Southern University, 1974 
M.S., Louisiana State University, 1978 


Submitted in partial fulfillment 
of the requirements for the degree of 


MASTER OF SCIENCE IN MECHANICAL ENGINEERING 
from the 
NAVAL POSTGRADUATE SCHOOL 
June 1996 





DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 


ABSTRACT 


The Structural Synthesis Transformation is used to conduct 
structural system identification in the frequency domain. For 
spatially complete cases where each of the frequency response 
functions at every degree of freedom of each of the coordinates of 
the modeled system are available it is shown that the theory exactly 
identifies all modeling errors. For spatially incomplete cases where 
the frequency response functions are available only at a proper 
subset of the degrees of freedom of the finite element model, single 
mode solutions are computed over intervals about the modes of the 
experimental system using matrices and complex valued line 
integrals. Methods of forming multiple mode solutions from the 


single mode solutions are explored. 
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I. INTRODUCTION 


The Finite Element (FE) method is a proven tool for modeling structural dynamic 
systems. As the complexity of the system increases the FE model may not accurately 
reflect the dynamic behavior of the system. To determine the extent to which the FE 
model accurately describes the physical system, a comparison of the dynamic response or 
modal parameters of the system as predicted by the FE model and the response or modal 
parameters of the physical system as determined by measurements of the dynamic 
response of the system can be made. Such a comparison can easily point out the 
differences in the dynamic behavior of the FE model and the physical system but fail to 
provide the necessary corrections to the FE model that will provide a more accurate FE 


model prediction of the dynamic behavior of the physical system. 


Structural system identification refers to procedures used to identify finite element 
modeling errors using dynamic test data. Localization is the process of identifying those 
degrees of freedom (DOF) of the FE model whose impedance differs from that of the 
physical system. We refer to this set of DOF as the error set. Identification is the process 
of finding matrices AK, AM, and AC that are corrections to the FE model stiffness, mass, 
and damping matrices. When corrections are installed, the frequency response of the 
corrected FE model better predicts the frequency response of the experimental system, 
and hence the modal parameters of the corrected FE model better predicts those of the 
experimental system. 


Structural system identification can be conducted using either modal or frequency 
domain methods. This thesis investigates a frequency domain method based on the 
structural synthesis transformation (SST) as outlined in Reference (1). When measured 
dynamic response data is available for each DOF of each coordinate of the FE model the 
identification is termed spatially complete and the SST can be used to exactly determine 
which elements of the FE model are in error. Additionally, the SST can be used to 
compute correction matrices AM, AK, and AC that can be used to correct the mass, 


stiffness, and damping matrices of the FE model. The dynamic response of the corrected 


FE model exactly matches that of the physical system. In the more frequent case that 
dynamic response data is not available at every DOF of the FE model, the identification is 
termed spatially incomplete. In this case the SST provides a frequency dependent 


solution. 


SST-based structural system identification uses the frequency response function 
(FRF) of the structure under consideration. The FRF is a quantitative description of the 
dynamic behavior of the structure. To obtain the FRF, known harmonic excitation forces 
are applied and the resulting harmonic response of the structure measured. The ratio of 
excitation forces to response at a coordinate evaluated at each frequency in a specified 


bandwidth defines the FRF. 


Hi. THEORY 


A. IMPEDANCE DESCRIPTION 
The impedance model of a given physical system can be defined by the relationship 


of structural displacement with respect to an applied force, 


tlle zl a 

The harmonic force and response vectors are denoted by "f" and "x" respectively. 

These vectors and the impedance matrix, Z, are in general complex-valued and frequency 
dependent. Subscripts "i" and "c" denote non-error and error DOF respectively. The 
superscript "a" denotes quantities calculated from a FE (analytic) model. If the values 


were obtained from experimental test data, the superscript would be "x". Thus, for the 


experimental model the impedance relationship would be given by: 


i Z mel, Nie 

i =| (2.1b) 
ip Lic Lee Xe 

In pratice, the elements of the experimental impedance matrix of Equation (2.1b) 


are unmeasurable quanties. To see this we expand row j of Equation (2.1b) to obtain 
Fj = ZX + Zj_Xq 00 +27, Xp (2.2) 


To measure an element z}, of Z’, requires that we impose a unit displacement at 


coordinate x; while physically holding all other other coordinates at zero displacements. 
This is not physically possible. Assuming, for purpose of definition, the availability of the 
experimental impedance matrix, the quantitative difference between the analytical and 
experimental systems, as a function of frequency, is described by the impedance error 


matrix. It is defined by the difference between the analytical and experimental impedance 


matrices. The error impedance matrix relationship is defined as: 


0 0 F y ZL ie 
= Zi Lic a wz ic (2 ; 3) 
OMINZ(O) lal Ze. 72 wide 72 


B. STRUCTURAL SYNTHESIS TRANSFORMATION 


Since the experimental impedance matrix, Z’, is in general unavailable, frequency 
domain structural synthesis is used to identify the impedance error matrix using FRF data 
exclusively. A structural synthesis transformation is constructed from AZ of Equation (2.3) 
which encompasses the FE model errors. This transformation is applied to the finite 


element model to produce an experimental system FRF. 

The FRF relates structural response to applied excitation. Given a FE model with 
impedance matrix, Z’, the FRF, H’, is the matrix inverse of the impedance matrix Z of 
equation (2.1a). For a static system ((QQ2=0) the FRF is simply the flexibility matrix (inverse 


stiffness matrix). Using the notation of Equation (2.1) we may partition H” as follows: 


x; Ay Ae hd: 
=| ji fi (2.4) 
x, ii ic el co if 
Generally, the "c" response coordinates experience applied forces due to both error 


impedances and externally applied forces, whereas "i" response coordinates experience 


only externally applied forces, such that, 


aire te (2.5a) 


and 


fat" (2.5b) 


Expanding Equation (2.4) and substituting the relations of Equation (2.5) yeilds 


the following relationships: 


p= Hf Pai eH, f* (2.6a) 
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x, = Hf" + Hef" +H fi” (2.6b) 


q 


If we include a copy of Equation (2.6b), in expanded matrix notation, Equation 


(2.6) reflects the three harmonic excitation terms to be considered, 1.e., 


x Ee, Heli 


x imma loo: alten ile. (2.7) 


Sli tet cell | 


Response coordinates "c" and "i", which are due to external forces, will hereafter 
be referred to as "e" coordinates, denoting their dependence on external force excitation. 


Consequently, the three excitation forces are condensed into two under the identity: 
(eh=[[47] | rl] (2.8a) 


{fF} = {07} (2.8b) 


Xe — ae ie Je 
bd Lae ae a as 


Equation (2.3) shows that the impedance error is defined as the difference between 


and Equation (2.7) reduces to 


the analytic and experimental impedance models. Hence, a transformation is required 


which uses the FRF relationship of Equation (2.9) to generate a similar relationship for the 


experimental system. The impedance error AZ provides the basis by which this 


transformation is developed. 


The impedance error matrix, AZ(QQ), is a function of frequency and satisfies: 
(7.} = -[az(ay.} 2.10 


where 


[AZ(Q)] = [[AK] + jQ[ AC] - 27[ AM] (2.11) 


for Q the forcing frequency, j = V/—1 and AK, AC, and AM stiffness, damping and mass 


error matrices comparable to those of the finite element formulation. The minus sign in 
Equation (2.10) reflects that the reaction forces imposed by impedance errors on the 


baseline model are being considered. Substituting the relationship 


Fale ee orem ey: 
Lis [o az. dct ans 
into Equation (2.9) yields: 
‘| {He Ale omit) _ 
b: TEA, We ; -az|te| 7% 
simplifying we get: 


eee (eS AZ (fala 
=| (2.14) 
X, H.. —H,,AZ |\x, 


cé 


Expanding Equation (2.14) into two equations and using "*" to denote a 


synthesized modified response results in 
x. = Hof, - H2AZx: (2.15a) 


eee AZ (2.15b) 


Rearranging Equation (2.15b) produces 
[I+ HzAZ|x; = His, 


Using the property of the frequency response function, 
x, = Heda 


[I+ Hz.AZ|x; = x, 
x. = [I+ HAZ] ‘x, 


Introducing Equation (2.16) into Equation (2.15b) results in 
x) = H2f,- HLAZ| I+ H2AZ] x, 


x = Hof, - H2AZ|[I1+ H2AZ) Hf 


Once again we recall the property of a FRF, 
x, = Hef, 


Combining Equations (2.17b) and (2.18a) yields 
H,, = He, - HZAZ|1+ HLAZ| He 


Noting that 
[1+ Heaz}" =[(az7 + H2)az] 
and applying the matrix property 
((a][5})° =[2) "fay" 
we get that 


H,, = He, - Hi[AZ" + He] He 


(2.16a) 


(2.16b) 


(2.16c) 


(2.16d) 


(2.17a) 


(2.17b) 


(2.18a) 


(2.18b) 


(2.19a) 


(2.19b) 


(2.20a) 


Replacing the superscript "*", which denotes the structures's synthesized coupled 


response, with the superscript "x" to indicate the test system response we arrive at 


Hz = Hz, - He[AZ" + Hz) He. (2.20b) 
In full matrix notation we have 
T 
H* H* H? #H? H? Pee 
ii ic | ii ieee ic [Az™ + Hé,| 1 ‘ (2.20c) 
Oe, pect laalelel a Tl ||aaleld ies 


Equation (2.20b) is the structural synthesis transformation equation. When the 
experimental system FRF, H” , is available the SST can be used to identify a frequency 
dependent impedance error matrix [AZ(Q)]. Additionally, using Equation (2.9), [AZ(Q)] 


can be decomposed into constitutent stiffness, mass, and damping errors. 


C. FREQUENCY DOMAIN LOCALIZATION 


We may rewrite Equation (2.20b) as 


AH. =e (2.21) 
where 
a Sea, (2.224) 
and 
D=[AZ" + Hi (2.22b) 


We define the localization matrix L as 
Eee. 2A SZ (2.23) 


using the expression of Equation (2.22a) in Equation (2.23) we can rewrite L as 


PZ eyez (2.24a) 


Expanding the "e" coordinate set into error and non error coordinates we get 


i ih Li Lic 
ze 7 


cc 


tl ees, asa" Su Zee 
[AP He Hel? ia (2.24b) 


Noting the the frequency response matrix is the inverse of the impedance matrix, L.e., 


OZ wt li. «ie I QO 
az ic fs ic = (2.25a) 
Ze de Bh kl. 0 / 
all mixed product coordinates of Equation (2.24b) must be zero and L simplifies to 


L= we Q=Qi 2.250 
|, i @ ~~] (2. ) 


D. ERROR IMPEDANCE 


We can solve Equation (2.20b) for the impedance error /4Z]. From Equation 
(2.25b) terms of Equation (2.20b) associated with non error coordinates may be assumed 


to be zero. We get the following form of the impedance error matrix 


[AZ] = ((#z]-[a2]) @ Q=Q; (2.26a) 
where 
[|= (ze) Lame] ) 2.268) 


Let S={Q), Qa, ... , Qa} be a set of frequencies where Q) < 02<...<Q,.1:<Qz,. If 
for each 1= 2,3, ..., n-1, we apply Equation (2.10) at each of the frequencies Qj), Q; , and 
Qi+1 and assemble the resulting three equations into a system of three equations in three 


unknowns, we get the matrix equation 


AZ(Q,,)| | 2 -Q2 jQ,.,1 1 AK, 


AZ(Q,) |=|I -Q7I  jQ,I | AM, (27%, 
AZ,(Q,41) I -0% I el AC, 


i+] 


Equation (2.27) can be used to decompose the frequency dependent impedance 
error into constitutent stiffness, mass, and damping error matrices, AK, AM, and AC at the 


frequencies Q;, i=2,3,....n-1. These constitutents matrices are in general frequency 


AK,(Q,) 
dependent. Denoting the solution of Equation (2.27) by , AM. (Q,) for i=2,3,...,.n-1 we 
AC,(Q,) 
obtain for each Q; i=2,3,....n-1, error stiffmess, mass, and damping matrices AK.(QQ), 
AMA(Q2;), and AC.(Q,) such that 


AZ,.(Q,) = AK,(Q,) — Q7AM,(Q,) + jQ,AC,(Q,) (2.28) 


The matrices AK.(92,), AM-(92;), and AC-(Q;) can be used to numerically correct 
the stiffness, mass, and damping matrices of the FE model at the frequencies Q); 
i=2,3,....n-1 in that the FRF of the corrected FE model at the frequencies Q; approximates 
the experimental system FRF at ©;. To express this symbolically, if.K’, Mé’, and C’ are the 
stiffness, mass, and damping matrices of the FE model and # is the FRF matrix of the 


experimental system at Q; 
H*(Q,) = (K? + AK.(Q,) -07(M? +AM,(Q,)) + j2,(C7 + AC.(Q,))) (2.29) 


where the "c" subscripted matrices are added at the corresponding error set coordinates of 
the "a" superscripted matrices. 
For a general set of frequencies, S={Q), Qy, ... , Qa}, we can form the system of n 


equations in three unknowns given by: 


AZ.(2)) I -QU ja 
; BRAK. 
42.(2,)|- I -Q?T jQ,1 | AM, 
: : KE 


c 


(2.30) 
| AZ, (2, )| , O72 {Mas 


AK. | | A4K.(5) 


Cc 


The solution | AM, |=| AM_(&)| of Equation (2.30) represents error stiffness, 


c 


AC. | | AC.(=) 


mass, and damping matrices which best corrects the FE model in a least squares sense. 


Equations (2.27) and (2.30) are fundamental to all that follows. 
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Ht. SPATIALLY COMPLETE STRUCTURAL IDENTIFICATION 


To illustrate the principles of SST based frequency domain structural identification 
we will use the FE model of a free-free beam. To simulate the experimental system we 
impose a 25% addition to the mass and stiffness of elements 3 and 4 ofa 10 element FE 
model. Figure 3-1 shows the finite element model of the beam and the spatially complete 
experimental system that results from imposing the mass and stiffness additions at element 
3 and 4 of the FE model. Table 3-1 shows the system frequencies of the analytic and 


experimental systems. 


oP ey PPh eg 


MASS OF ELEMENT 3 = 1.954X10*-6 LBM 
MASS OF ELEMENT 4 = 1.954X10*-6 LBM 
STIFFNESS OF ELEMENT 3 = 417.88 LBFAN 
STIFFNESS OF ELEMENT 4 = 417.88 LBF/AN 


2 tt 
REE RETEES:E 


ANALYTIC SYSTEM 

NUMBER OF ELEMENTS = 10 

ferent e co cos faa) ROTATIONAL DOF AT NODE 

El = 74500 LBF-IN*2 

SECTIONAL AREA = 0.8359 IN“2 - TRANSLATIONAL DOF AT NODE 
WEIGHT DENSITY = 0.05 LBFAN*3 

ELEMENTAL MASS = 1.563X10*-6 LBM 

ELEMENTAL STIFFNESS = 334.3 LBF/AN 





Figure 3- 1 Spatially complete Analytic and Experimental Systems. 
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MODE ANALYTIC (Hz) EXPERIMENTAL 


p  54E 41 P2442 42 





Table 3-1 System Frequencies of spatially complete Analytic and Experimental systems. 


We will denote by /", K°, and C* the mass stiffness and damping matrices of 
the FE model and by M@*, K*, and C* the mass stiffness and damping matrices of the 


experimental system. The impedance matrices of the analytic and simulated experimental 


systems are given by: 
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Z*(Q) = K* + jQC* -Q?M’ (3.1a) 
Z*(Q) = K* + jac* -Q’°M* (3.1b) 


The FRF matrices of the analytic and simulated experimental systems are given by: 


HX@) = 740) ' (3.2a) 
H*(Q) = 27(Q)" (3.2b) 


Figure 3-2 shows a comparison of the driving point FRF at DOF 1 of the analytic 


and experimental system. 
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Figure 3-2 Analytic FRF vs simulated Experimental FRF. 


14 


800 


1000 





1200 





For a given frequency Qo, using Equation (2.24a), we can form the localization 
matrix L. Plotting the diagonal elements of L versus the associated DOF we obtain the 
plot shown in Figure 3.3. For each frequency, Q, in a given frequency range, we can 
compute the diagonal of L at © . Assembling the diagonals over the frequency range of 
our system into a rectangular matrix and performing a MATLAB mesh plot of the 
resulting rectangular matrix we obtain the surface plot shown in figure 3.4. Figure 3.4 
shows the frequency dependency of the localization matrix diagonals. As our system is 
spatially complete Equation (2.25) forces all non error coordinates to be zero. From 
Figure 3-3 we can determine that the locations of the nonzero diagonal values are DOF 9, 


10, 11, and 12. We denote by C, the set of DOF for which L()) is non zero. For our 
beam system C,, = {9,10,11,12} . Figure 3.5 shows the frequency dependence of typical 
error and non error set diagonal elements of the localization matrix L. 

Using the set, Cex, which results from the localization, we can perform a 
partitioning of the FRF matrix as described by Equation(2.4). We can now apply 
Equation (2.26) to compute AZ as a matrix function of frequency over the frequency 
range of our system. We then use Equation (2.27) to decompose AZ into its constituent 
components AK, AM, and AC. We get exact solutions of error stiffness, mass, and 


damping as shown in figures 3-7, 3-8 and 3-9. The MATLAB Routine SST.M of 
Appendix A can be used to accomplish the above steps. 


For each Q in the frequency range of our system we can form the sum: 
Z corr (2) = Z7,(Q) + AZ(Q) (3.3) 
We refer to 
Hoon = (Zeor) (3.4) 


as the corrected FRF of the analytic model. Ho, is the FRF of the model that results from 
installing the corrections as identified by the Equation (2.26a). Figure 3-10 shows a 
comparison plot of the FRF of the corrected model and the FRF of the experimental 


system. Figure 3-10 clearly shows the exactness of the SST solution in the case of a 


1S 


spatially complete system. The experimental and corrected model FRFs are identical to 
within plot resolution. 


At each frequency, Q, in the frequency range of our system we can form the sum 


Zo" (Q) = Z2.(Q) + AK(Q) + fAAC(Q) - Q?AM(Q) (3.5) 
We refer to: 
Hepa) (3.6) 


as the constituent corrected FRF. Figure 3-11 is a comparison plot of the constitutent 
corrected FRF and experimental for our system. The offset between the two plots is 


exactly equal to the sampling frequency AQ. 
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Figure 3- 3 Spatially complete localization matrix diagonal at Q=352 Hz. 
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Figure 3- 4 Frequency dependence of spatially complete localization matrix diagonals. 
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Figure 3- 5 Frequency dependence of spatially complete localization matrix error set 


DOF. Units are Ibf/in (Log10 of). 
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Figure 3- 6 Spatially complete Analytic Impedance vs Experimental impedance. 
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Figure 3- 7 Computed stiffness vs true stiffness for spatially complete beam. Plots are 


identical within plot resolution. 
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Figure 3- 8 Computed mass vs true mass for spatially complete beam. Plots are identical 


within plot resolution. 
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Figure 3- 9 Computed damping vs true damping for spatially complete beam. Plots are 


identical within plot resolution. 
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Figure 3- 10 Spatially complete Experimental FRF vs AZ corrected FRF. Plots are 


identical within plot resolution. 
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Figure 3- 11 Spatially complete experimental FRF vs AK, AM, and AC corrected FRF. 
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IV. SPATIALLY INCOMPLETE STRUCTURAL 
IDENTIFICATION 


A. GENERAL DESCRIPTION 


To illustrate SST based frequency domain structural identification when the 
physical system under consideration is spatially incomplete we will again use the FE model 
of a free-free beam. We simulate the experimental system by imposing a 25% addition to 
the mass and stiffness of elements 3 and 4 of a 10 element FE model. Figure 4-1 shows the 
FE modeled beam and the spatially incomplete system that results if FRF data is available 
only at the displacement DOF of the FE system. 


EXPERIMENTAL SYSTEM 

MASS OF ELEMENT 3 = 1.954X10-6 LBM 
MASS OF ELEMENT 4 = 1.954X10*-6 LBM 
STIFFNESS OF ELEMENT 3 = 417.88 LBFAN 
STIFFNESS OF ELEMENT 4 = 417.88 LBFAN 


NUMBER OF ELEMENTS = 10 
LENGTH = 60.625 

El = 74500 LBF-IN*2 q 

SECTIONAL AREA = 0.8359 IN*2 TRANSLATIONAL DOF AT NODE 
WEIGHT DENSITY = 0.05 LBFAN“3 

ELEMENTAL MASS = 1.563X104-6 LBM 

ELEMENTAL STIFFNESS = 334.3 LBF/AN 


f ROTATIONAL DOF AT NODE 





Figure 4- 1 Analytic and experimental spatially incomplete systems. 
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To obtain a simulated FRF for our spatially incomplete beam we denote by 
M°’, K*, and C°’ the mass stiffness and damping matrices of the FE model and by M“, 
K*, and C” the mass stiffness and damping matrices of the experimental system. M~, 


K*,and C”* are obtained by imposing 15% mass and stiffness and a 15% mass additions 
to the elemental matrices of element 5 and 6 respectively of the FE model. The impedance 
matrix of the simulated spatially complete beam is given by 

Z*(Q) = K* + jQC* -Q°M (4-1) 
The FRF matrix of the simulated spatially complete system is given by 


ao) = 7G (4-2) 


We introduce the terminology Analysis set and Omitted set where the Analysis set 
(A-set) is that set of DOF for which experimental FRF data is available and the Omitted 
set (O-set) is that set of DOF of the experimental system for which experimental FRF data 
is unavailable. For our simulated experimental system the A-set consist of the odd 
numbered translational DOF and the O-set consists of the even numbered rotational DOF, 


1.€., 


A —ser= {1,3,5;...,21; (4.3a) 
O- set = {2,4,6,...,22} (4.3b) 


We obtain the simulated FRF of our spatially incomplete beam by physically 
extracting the rows and columns of the simulated spatially complete matrix, H’, for which 
FRF data would be available. In our system these are the rotational DOF and all even 
numbered rows and columns are omitted from the simulated spatially complete FRF to 


obtain a spatially incomplete FRF which we will denote by H". 
For a fixed Qo, H"(Q)) is a square matrix of size (length(A-ser)) by (length(4- 
set)). For our FE model as currently defined, H’, is of size (number of DOF) by (number 


of DOF) and is, as is true of most real world cases, of larger size than H" . In order to 


employ the structural synthesis transformation we need to reduce the size of H” to that of 


74g) 


H’ . To this end we will consider two reduction methods, FRF matrix extraction (Ref. 1], 
and the Improved Reduced System as given in [Ref. 2]. 


B. EXTRACTION REDUCTION METHOD 


In order to reduce H” by the extraction method we simply extract from the full 
order H those rows and column which correspond to A-set coordinates. Partitioning the 


impedance and full FRF matrices of our analytical system 


Ze 
Z; -| ~ 4 (4.4a) 
eer 2, 
ae 
Feel. (4.4b) 
LT cc" ldeo 
and using the identity 
Loe Le Wie Tl I 0 
CE ti Ali oe (4.5) 
a LW, Td, Or Tf 
we obtain the relationship 
Hyg = (Zan - ZeoZos Zea) (4.6) 


We denote the reduced Analytic FRF of Equation (4.6) by H’. 


C. O-SET SYSTEM 


Equation (4.6) relates the extracted reduced order FRF of the analytic model to the 
impedance of the full order order model. Taking advantage of the identity 


[Z,.] | = (det[Z,.]) adi[Z.0] (4.7) 


and replacing Zo. by Kootj2Coo-j§2%Moo (where det(e) and adj(e) represent the 
determinant and adjount respectively), we see that an element H?(Q) is large for those 


frequencies 92, where £2, is an eigenvalue of the O-set system, the O-set system being that 
FE model having stiffness, mass and damping matrices Ky., Moo, and Cy. respectively. 
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D. IMPROVED REDUCTION SYSTEM 


To use the improved reduction method (IRS) we first partition Z’ using the A and 
O sets, then adjust Equation (2.1) to reflect this new coordinate system obtaining 


Ja| _| Zee Zao |) *e (4.8) 
5 calala rena lcoaldies 
Expanding Equation (4.8) ito two equations yields 
Wo— 2. xX, + Zio (4.9a) 
i= Ze (4.9b) 


O-set coordinates are not associated with FRF data measurement locations on the 
physical structure, therefore the forcing function at O-set coordinates can be set to zero. 
Making these substitution into Equation (4.9b) and solving for the generalized structural 


response coordinates leads to: 


x, =-Z1Z%x (4.10a) 


00°" Oa" a 


Xq| _ I 
} = _ atte (4.10b) 


Substituting these results into Equation (4.8) yields, 


Ame, Zz I 
WE S[k a 


hence 
{fa} =[Zoe - ZaoZoo Zea |{Xa } (4.12) 


When 02=0 the Equation (4.10a) yields the static reduction relationship between 
omitted and retained coordinates and is given by 


{xo} =[-Kio Koal {xa} (4.13) 


Zo 


The IRS relationship 1s given by 
{Xo} =|—Kap Koa + TMgaeK sam |{%a} (4.14) 


where 


tT Sieie= KM KeLK 


00 oa 


(4.15) 


and Kgat and Mga are the statically reduced [Ref. 2, 3] stiffness and mass matrices. 

Unlike the spatially complete case where we only had to consider two system (the 
analytic and the experimental systems), in the case of a spatially incomplete system there 
are actually five systems with which we must concern ourselves [Ref 4]: the analytic 
system, the experimental system, the reduced analytic system that results from conducting 
dynamic reduction on the mass and stiffness matrices of the analytic system, and the 
omitted systems of both the analytic and the experimental system. Table 4-1 shows the 


frequencies of each the first four of these systems. 


Figure 4-2 shows a comparision of the analytic and experimental FRF of our 
spatially incomplete beam. We see that by using IRS reduction of the analytic system those 
modes above approximately 1000 Hz i.e., those modes associated with the reduced out 
rotations, are not present. Figure 4-3 shows a similar comparision where we have used 
extraction reduction and the higher modes are present. In all that follows we shall use IRS 
reduction of our analytical systems. Figure 4-4 shows the localization matrix diagonal at 
O=196.1 Hz while Figure 4-5 shows the frequency dependency of the localization 
diagonals over the frequency range our our spatially incomplete beam. For a spatially 
incomplete system the determination of the locations of the error coordinates is not a 
clearly defined task. We shall not discuss the problem of actually determining the exact 
error set in the spatially incomplete case and will use our knowledge of the true location of 
the error coordinates to aid in our localization. As we shall make use of the concept of the 
size of the error set again, we will simply note that using a reduction method like IRS 
causes errors to be 'smeared' in the reduced analytic model. Table 4-2 shows the A-set and 


O-set of our simulated beam along with the locations of the true errors as well as the 
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computed C-set at Q=196 Hz. The computed C-set is the set of all DOF having a diagonal 


entry whose absolute value exceeds a given tolerance. 


Mode Exp (Hz Reduced | O set (Hz 
(Hz) 


res ar 
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Table 4-1 Analytic, Experimental, Reduced, and Omitted System Frequencies of Spatially 






















Incomplete Beam. 
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Table 4-2 Analytic set, Omitted set, Computed C set, and True C set DOF for spatially 






incomplete beam. 


32 





H(7,7) in/lbf (Log10 of) 


0 200 400 600 800 1000 
Omega (Hz) 


Figure 4- 2 Analytic FRF vs experimental FRF for spatially incomplete beam using IRS 


reduction. 
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Figure 4- 3 Analytic FRF vs experimental FRF for spatially incomplete beam using 


extraction reduction. 
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Figure 4- 4 Localization matrix diagonal at Q=196.1 Hz for spatially incomplete beam. 
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Figure 4- 5 Frequency dependence of localization matrix diagonals for spatially 


incomplete beam. 
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Figure 4- 6 Frequency dependence of error and non error set localization matrix 


diagonals for spatially incomplete beam. Units are lbf/in (Log 10 of). 
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Figure 4- 7 Analytic vs experimental impedance for spatially incomplete beam. 
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Figure 4- 8 Computed Stiffness vs True Stiffness for spatially incomplete beam. 
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Figure 4- 9 Computed Mass vs True Mass for spatially incomplete beam. 


40 





DC(7,7) lbf-sec/in (Log10 of) 


0 200 400 600 800 1000 
Omega (Hz) 


Figure 4- 10 Computed Damping vs True Damping for spatially incomplete beam. 
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Figure 4- 11 AZ Corrected FRF vs experimental FRF for spatially incomplete beam. Plots 


are identical to within plot resolution. 
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Figure 4- 12 AK/AM/AC Corrected FRF vs experimental FRF for spatially incomplete 
beam. The offset of the two plots is equal to the sampling frequency AQ. 
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V. SINGLE MODE SOLUTIONS 


A. SINGLE MODE MATRIX SOLUTIONS 


In the case of a spatially complete system we have seen that the SST yielded 
frequency independent error matrices AK, AM, and AC that could be used to correct the 
stiffness, mass, and damping matices of the FE model in such a manner that the FRF of the 
corrected FE model was exactly equal to that of the experimental system. In the case ofa 
spatially incomplete system the constituent solutions AK, AM, and AC given in chapters 
III and IV were in general frequency dependent solutions which serve only as corrections 
to the reduced FE model. Ideally we seek frequency independent solutions which are 
corrections to the full FE model. For now we shall only deal with the simpler problem of 
trying to find frequency independent solutions which serve as corrections to the reduced 


FE model. 


We shall employ Equation (2.30). For a given experimental system mode, @,, 
consider a frequency bandwidth [Q), Quy] such that @,€[Q), Qu]. Let =={OQ) Op,..., Qn} be 
a frequency sampling of the bandwidth [Q), Quy], 1.e., Q:S$Q;<Q, for each QieS. For each 
Qie= we can form the error impedance matrix AZ(Q,). We apply Equation (2.30) to the 
partition yeilding the following system of m equations in three unknowns: 

AZ(Q,)| [I -Q?21 jQl 
ne 
AZ,(Q,) |=|2 -Q?2 jQ,1 | AM® (5.1) 
AGE 
AZ,(2,,)4) |¥ -Q21 jor 

We will demonstrate by example that for properly chosen bandwidths [Q), Qu], the 
solution AK, AM’, AC?" of Equation (5.1) approximately corrects the Analytic FRF 
over the frequency bandwidth [Q),Q,], i.e., if H’(Q) is the value of the FRF matrix of the 


experimental system at a frequency Q, H(Q) the value of the FE reduced analytic FRF 
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matrix at Q, and H%,,(Q) the value of the single mode corrected reduced Analytic FRF 


corr 


matrix at Q for an experimental system mode @, where Qe [12),Q.] then 





H"(Q) - H,,(Q)|s|H"(Q) - #°(Q) (5.2) 


We will refer to the solution AK, AM?", ACY" as a single mode solution at @n. 


In the case of a spatially complete system, such as the spatially complete beam discussed in 
chapter III, all single mode solutions are found to be identical to the unique frequency 
independent solution that was obtained in chapter III, hence the corrections are valid 
throughout the frequency range of the experimental system and the left hand side of 


Equation 5.2 is zero. For a spatially incomplete system this is not the case. 


In support of Equation (5.2) we chose a 1 Hz frequency bandwidth centered on 
mode 1 (@,=25.21 Hz) of our spatially incomplete beam as defined in Chapter II. For our 
frequency sampling, =, we will use a sampling frequency of .5 Hz which results in 3 points 
which are equally spaced over the interval, =={24.7178 Hz, 25.2178 Hz, 25.7178 Hz}. 
Solving Equation (5.1) we get 

—3.91 -13.60 21.70 


AK =|-13.60 170.5 -187.5 (5.3a) 
-21.70 -187.5 194.6 


—0.0001 0.0006 -—0.0005 
AM* =| 0.0006 -0.0033 0.0026 (5.3b) 
—0.0005 -0.0026 -0.0016 


0+0.02707 0-0.2291j 0+0.2277; 
AC” =|0-0.2291j 0+1.50787 01.4627; (5.3c) 
0+0.22777 0-146277 0+13781) 


Figure 5-1 shows a comparison of experimental, uncorrected and mode 1 
corrected FRFs for our spatially incomplete beam using a 1 Hz frequency bandwidth 
centered on Mode | and a frequency sampling consisting of 3 frequencies equally spaced 


over the bandwidth. Figure 5-2 shows the results of including mode 1 and mode 2 in the 


4§ 


frequency bandwidth [Q), Qu]. As Figure 5-2 shows Equation 5.1 is not as accurate when 
multiple frequencies are included in the bandwidth. Figure 5-3 is a comparison plot of the 
experimental and analytic FRF versus the single mode matrix solution corrected FRFs over 
25, 10, and 1 Hz bandwidths using a 3 point frequency sampling. Figure 5-3 shows that 
the precedure is sensitive to the size of the bandwidth over which the solution is 
computed, i.e., better accuracy is achieved with smaller bandwidths. Figure 5-4 is a 
comparison plot of the experimental and uncorrected analytic FRFs versus the single mode 
matrix solution corrected FRFs computed over a 1 Hz bandwidth using frequency 
samplings of 200, 50, 10, and 3 points equally spaced over the bandwidth. Figure 5-4 


shows that the procedure is fairly insensitive to sampling frequency. 


Figure 5-5 is a comparison plot of the experimental and analytic FRFs versus the 
single mode solution corrected FRF for a 1 Hz bandwidth using a 3 point frequency 
sampling in the case of a spatially complete beam. As Figure 5-5 shows the procedure 1s 


exact for spatially complete systems. 
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Figure 5- 1 Experimental and uncorrected Analytic FRFs vs single mode solution at 
mode 1| corrected FRF using a 3 point frequency sampling of a 1 Hz bandwidth. 
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Figure 5- 2 Experimental and uncorrected analytic FRFs vs single mode corrected FRF 
using a 15 point frequency sampling of a bandwidth that includes modes 1 and 2. 
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Figure 5- 3 Spatially incomplete experimental FRF vs single mode matrix solutions at 


mode 1 using 3 point frequency samplings of 25, 10, and 1 Hz bandwidths. 
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Figure 5- 4 Spatially incomplete experimental FRF vs single mode matrix solutions at 


mode 1 using 3, 10, 50, and 200 point frequency samplings of a 1 Hz bandwidth. 
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Figure 5- 5 Spatially complete experimental FRF vs single mode matrix solutions at 


mode 1 using a 3 point frequency samplings of a 1 Hz bandwidth. 
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B. SINGLE MODE INTEGRAL SOLUTIONS 


In what follows we wish to employ the integral formulas of reference (5) to obtain 
a results similar to that of Equation (5.2). To accomplish this we shall need to recast the 


impedance equations of Chapter II in terms of velocity. We start with the equation of 


motion of a general 2nd order linear system 
Kx +Cx+ Me = f 


where 


x = Xe™ 
f= Fe™ 


If we differentiate Equation (5.4b) we get that 
X = jAXe™ = jOx 


hence 


differentiating Equation (5.4b) twice we get that 
¥ = jA(jOXe™) = jQx 
substituting Equations (5.6 )and (5.7) into Equation (5.4a) we get 
K Ee + Cx+ MjQx = f 
JQ 
We can rewrite equation (5.8) as 


Es c+} =a 
JQ 


(5.4a) 


(5.4b) 


(5.4c) 


(5.5) 


(5.6) 


(5.7) 


(5.8) 


(5.9) 


Following reference (5) we can take the Laplace transform of equation (5.9) and 


thus write the impedance of the linear system in terms of the complex Laplace parameter s 


as 


S2 





2(s) = Ms+C+= (5-10) 


In reference (5) the impedance matrix of a general system is represented as an 
infinite Laurent series about the orgin in powers of the complex Laplace parameter s as 


Z(s) = > 4,5" = Lim,_,.(4,8 +°:+A,S + Ay a atte nt ii as (5.11) 
5S S 


n=-—@ 


After defining a truncated Laurent series 


Ta As OA nee eee (5.12) 
S S 


which approximates the impedance matrix, an error function £(s), and a cost J are defined 


as 


E(s) = Z(s)- Z(s) (5.13) 
and 


J=q_|A(s)||,\7(s)- as (5.14) 


where W(s) is a complex valued weighting function, the subscript E denotes the euclidian 
norm of a matrix and the integration is performed over a prescribed path P in the complex 
plane. By setting the partial derivatives of J with respect to M,K,andC to zero the 
authors of reference (5) obtained expressions for matricesM,K, andC which 
approximate the stiffness, mass and damping matrices about a mode of the system. The 


expression for M,K,and C areas follows: 





— 1 | ‘hk: 
M= a agp hoz, (iQ)|W(Q)|dQ - ‘he Z,(iQ)|W(iQ)|dQ 
(5-15a) 
= _ 1 . 
Gre 1 fz, (aypr(njan (5-15b) 


2, 
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“fae... lee 
K=—— “| QZ, (iQ)|W(iQ)|da - 2) 5% (iQ)|W(iQ)|dOo 
(5-15c) 
where 
Q, 
a= [Q?|W(iQ)dO (5-15d) 
O, 
“2 I 
b= ! orl (iQ) ao (5-15e) 
2, 
c= ||W(i2)"do (5-15f) 
Z(iQ) = Z,(iQ) + iZ,(iQ) (5-15g) 


We will apply Equation (5.15) to AZ as defined by the SST and obtam 
matrices AK , AM, and AC which will serve as correction matrices of the FRF of the 
system about a given mode of the system. 

As an example, let us use Equation (5.15) to compute the (1,1) element of the 
correction matrix AC at mode 1 (@;=25.21 Hz) for our spatially incomplete system. We 
take the weighting function to be Wiif2)=1/7ig2) , the path P 1s taken to be along the 
imaginary axis from 155.30 rads/sec to 161.58 rads/sec. 

To compute this line integral we shall use the trapezoid rule with a three pomt 
frequency sampling 

= = {155.30rads / sec,158.44rads / sec,161.58rads / sec} (5.16) 
of the straight line path P (sampling frequency of x rads/sec). We first compute the value 


of the constant a of Equation (5.15d). Using vector notation we have that 
Q(&) = [155.30rads / sec,158.44rads / sec,161.58rads / sec] (5.17a) 


and 


W(S) = [0 — 0.0064 j sec/ rads,0 — 0.0063) sec/ rads,0 — 0.0062 sec/ rads} (5.17b) 
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hence the integrand is the vector 
2? (E) -|W(S)| = [155.30,158.44,161.58] (5.18) 


Using the trapezoid rule to compute the integral, we have that 


161.58 
ieee (155.30 161.58 
a= { 2°(5)|W(B|d2- 7 +158.44 4 aoe = 995,56 rads he 2 (5.19) 


oe. op Z 
b=.00000158 sec’/rads’ and c=0.0397 are computed in a similar manner. 

To actually compute the matrix AC ,we use the SST to compute the matrices 
AZ(155.30 rads/sec), AZ(155.44 rads/sec), and AZ(161.58 rads/sec) for the points of the 
frequency sampling of the straight line path P. Using Equation (2.26) the matrix 
AZ(155.30 rads/sec) 1s seen to be 


(0.64E —9)+0.02j (-0.04E-8)-0.04j (0.03E-7)+0.01) 


Az(155.44rads/_ |) =| (-0.4E-8)-0.04j (-04E-9)-0.04j (0.25E-7)-010; |"I/ 


(03E-8)+0.01j (0.25E-7)-0.10j (024E-7)+015; 
(5.20) 


The integral in Equation (5.15b) 1s computed on an element by element basis. For 
the (1,1) element of Equation (5.15b), if we collect the (1,1) elements of the matrices 
AZ(155.30 rads/sec), AZ(155.44 rads/sec), and AZ(161.58 rads/sec) and use the 
weighting vector given in Equation (5.17b) the integrand of Equation (5.15b) is the vector 


AZ (=)|W(S)| = [0.0064 0.040 -0.848]0.64 0.63 0.62(E-9)L/ sec/ 








(5.21a) 

AZ (S)|W(S)| = [0.0041 02885 0.0036(E-9)/2// sec/_ (5.21b) 
Computing the integral we have 

AC (1,1) = a AZ) (5) W(S)ldQ = star z)a + 0.2885 + oO \CE ao) 

(5.22a) 

AC (1,1) = 0.02328 - 6 / (5.22b) 


Performing the above procedures for the remaining elements of Equation (5.15b) we get 
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l 


n 








oe eo ——> <a 


=> > Ne 


lc —_— = = 





> 4+ = = = 





S i> 
_ > = 
= 








0.0232 0.0443 0.0300 
AC =|-0.0443 0.1278 -0.0914(E-6) YY (5.23a) 
0.0300 -0.0914 0.0645 


In a similar manner we have 


0.0001 -0.0001 0.0002 
AM =|-0.0001 0.0015 0.0020 \/bm (5.23b) 
0.0002 -0.0020 0.0027 


—6.0476 4.5347 3.6757 
AK =| 4.5347 51.2196 -71.7061|'"7/, (5.23¢) 

3.6757 -71.7061 85.4417 
Figure 5-6 is a comparison plot of experimental and uncorrected analytic FRFs 
versus single mode integral solutions at mode 1 using a frequency sampling of 3 points 
equally spaced over bandwidths of 25, 10, and 1 Hz. As with the matrix solutions, the 
integral solutions are sensitive to the length of the bandwidth over which the solutions are 
computed with smaller bandwidths yielding more accurate solutions. Figure 5-7 is a 
comparison plot of experimental and uncorrected analytic FRFs versus single mode 
integral solutions at mode 1 over a 1 Hz bandwidth using frequency samplings of 3,10,50, 
and 200 points from the bandwidth. As with the matrix solutions, the integral solutions are 


also insensitive to the sampling frequency. 


Figures 5-8, 5-9, and 5-10 are comparison plots of matrix and integral solutions at 
mode 1 computed over 25, 10, and 1 Hz bandwidths respectively, using three point 
frequency samplings. These three plots show that for these frequency sampling, the matrix 
solutions are more accurate than the integral solutions. We must note that we have used 
the trapezoid rule for computation of the integrals. It is expected that better accuracy 
would be achieved from the integral solutions if a higher order integration method such as 


Simpson's rule were used. 
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Figure 5- 6 Spatially incomplete experimental FRF vs single mode integral solutions at 


mode 1 using 3 point frequency samplings of 25, 10, and 1 Hz bandwidths. 
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Figure 5-7 Spatially incomplete experimental FRF vs single mode integral solutions at 


mode 1 using 3, 10, 50, and 200 point frequency samplings of a 1 Hz bandwidth. 
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Figure 5- 8 Experimental FRF vs single mode integral and matrix solutions at mode 1 


using a 25 Hz bandwidth with a 3 point frequency sampling. 
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Figure 5- 9 Experimental FRF vs single mode integral and matrix solutions at mode 1 


using a 10 Hz bandwidth with a 3 point frequency sampling. 
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Figure 5- 10 Experimental FRF vs single mode integral and matrix solutions at mode 1 


using a 1 Hz bandwidth with a 3 point frequency sampling. 
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VI. MULTIPLE MODE SOLUTIONS 


A. MULTIPLE MODE MATRIX SOLUTIONS 


We wish to extend the results of Chapter V to multiple modes of the system under 
consideration. To accomplish this we shall simply extend Equation (5.1) over multiple 
frequency samplings. For a set of system modes {@), @,...,@:}, let E(@)={Qir, 
Qi2,..-,Qix} be a frequency sampling of a bandwidth [Q)),Qi,] about ; for 1=1,2,....n. We 
shall apply Equation (2.30) to the frequency sampling 


==U,..5lo,) (6.1) 


which is the set theoretic union of the sets =(@j;) for =1,2,...,.n. If, for simplicy, we restrict 
ourselves to equally sized frequency samplings we obtain a set on nk equation in three 
unknowns 

AZ,(Q,;) I =oi/ JX L 


est ; 6 AKS 
AZAQ| | =| tie Otay Oe AM; (6.2) 
. AC? 
AZ,(2,) LI On JQ L | 
AK? 
One can solve Equation (6.2) for = | as we have done in Chapter V using the 
AC= 


familiar MATLAB puedoinverse function but better results are obtained if we weight the 
equations as discussed in [Ref. 6]. A good weighting is to assign the weight @; to those 
equations associated with the frequency sampling =(«a;). Alternately we could assign the 
weight 1/a@; to those equations associated with the frequency sampling =(@;). The @; 
weighting results in the solution being more accurate at the lower modes while the 1/; 


weighting gives better accuracy at the higher modes. 


We note that the solution to Equation (6.2) approximately corrects the analytic 
FRF in the sense of Equation (5.1) when the matrices AZ(Q,,) as given by Equation (2.26), 
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are computed over a ‘good’ error set, Cex. In general, due to the smearing that occurs 
during reduction of the full analytic system to the reduced analytic system, the set Ce, is 
not the same set as the intersection of the spatially complete error set and the A set. It is in 
general a larger set than this intersection and the size is a function of the type of error, i.e., 
mass, damping, or stiffness error and the locations of the errors. Figure 6-1 is a 
comparison plot of experimental and uncorrected analytic FRFs versus the multiple mode 
matrix solution corrected FRF with the solution having been computed over modes | and 
2 using a 1 Hz bandwidth about each mode with a 3 point frequency sampling over the 
bandwidth at each of the modes. Figure 6-2 is a comparison plot using a matrix solution 
computed over mode | through 4 again using a 1 Hz bandwidth about each mode with a 3 


point frequency sampling over the bandwidth at each of the modes. 
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Figure 6- 1 Experimental FRF vs multiple mode matrix solutions at modes 1 and 2 using 


a 1 Hz bandwidth with a 3 point frequency samplings at each mode. 
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Figure 6- 2 Experimental FRF vs multiple mode matrix solutions at modes 1 through 4 
using a 1 Hz bandwidth with a 3 point frequency samplings at each mode. 
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B. MULTIPLE MODE INTEGRAL SOLUTIONS 


We wish to extend single mode integral solutions to multiple modes of the system. 
As with the multiple mode matrix solutions we simply take as the path of integration a 
path which is the set theoretic union of suitable paths at each of the modes under 
consideration. As with the matrix solutions one can chose the weighting function W(s) to 
be 1/(Q) or Q to achieve lower or higher mode accuracy. As stated in the previous section 
it is required that a 'good' error set is known. Figure 6-3 is a comparison of experimental 
and uncorrected FRF versus multiple mode integral solution corrected FRF with the 
solution having been computed over modes 1 and 2 using a 1 Hz bandwidth at each mode 
with a 3 point frequency sampling over the bandwidth at each of the modes. Figure 6-4 is 
a comparison using a solution computed over mode 1 through 4 again using a 1 Hz 
bandwidth at each mode with a 3 point frequency sampling over the bandwidth at each of 
the modes. Figure 6-5 is a comparison plot of the multiple mode matrix and integral 
solutions over modes 1 through mode 4 using a 1 Hz bandwidth at each mode with a 3 


point frequency sampling over the bandwidth at each of the modes. 
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Figure 6- 3 Experimental FRF vs multiple mode integral solutions at modes 1 and 2 using 
a 1 Hz bandwidth with a 3 point frequency samplings at each mode. 
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Figure 6- 4 Experimental FRF vs multiple mode integral solution at modes | through 4 
using a 1 Hz bandwidth with a 3 point frequency samplings at each mode. 
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Figure 6- 5 Experimental FRF vs multiple mode matrix and integral solutions at modes 1 


and 2 using a 1 Hz bandwidth with a 3 point frequency samplings at each mode. 


69 





C. SINGLE POINT MULTIPLE MODE SOLUTIONS 


The results of sections A and B can be performed using partitions consisting of a 
single point. Figure 6-6 shows the results of performing a multiple mode matrix solution 
over the first three modes of our spatially incomplete beam using a single point partition at 
each of the modes under consideration. Figure 6-7 shows the results of performing a 
multiple point integral solution over the first three modes of our spatially incomplete beam 
using a single poit partition at each of the modes under consideration. Figure 6-8 
compares the multiple mode matrix and integral solution computed over the first three 
modes of our spatially incomplete beam with the analytic FRF. In figure 6-9 we compare a 
multiple mode matrix solution computed over the first three modes of a spatially complete 
beam using a single point partition at each of the modes with the beam's spatially complete 
analytic FRF. In figure 6-10 we compare a multiple mode integral solution computed over 
the first three modes of a spatially complete beam using a single point partition at each of 
the modes with the beam's spatially complete analytic FRF. Figure 6-11 compares multiple 
mode matrix and intgral solution computed over the first three modes of a spatially 
complete beam with the beam's spatially complete analytic FRF. Figures 6-9 and 6-10 


show that the single point multiple mode solutions are exact. 
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Figure 6- 6 Experimental FRF vs multiple mode matrix solution computed at modes 1 
through 3 using a 1 Hz bandwidth with a single point frequency samplings at each mode 


for a spatially incomplete beam. 
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Figure 6- 7 Experimental FRF vs multiple mode integral solution computed at modes 1 
through 3 using a 1 Hz bandwidth with a single point frequency samplings at each mode 


for a spatially incomplete beam. 


72 





H(9,9) in/lbf (log10 of) 


——— Experimental FRF 
‘|— — Corrected INT Anal FRF 
— — Corrected MAT Anal FRF 
-—--— Uncorrected Anal FRF 
tr Included Modes 





0 20 40 60 80 100 120 140 160 180 
Omega (Hz) 


Figure 6- 8 Experimental FRF vs multiple mode matrix and integral solutions computed at 
modes 1 through 3 using a 1 Hz bandwidth with single point frequency samplings at each 


mode for a spatially incomplete beam. 
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Figure 6- 9 Experimental FRF vs multiple mode matrix solution computed at modes 1 
through 3 using a 1 Hz bandwidth with a single point frequency samplings at each mode 


for a spatially complete beam. 
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Figure 6- 10 Experimental FRF vs multiple mode integral solution computed at modes 1 
through 3 using a 1 Hz bandwidth with a single point frequency samplings at each mode 
for a spatially complete beam. 
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Figure 6- 11 Experimental FRF vs multiple mode matrix and integral solution computed 
at modes 1 through 3 using a 1 Hz bandwidth with a single point frequency samplings at 


each mode for a spatially complete beam. 
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Figure 6- 12 Experimental FRF vs multiple mode matrix solutions computed over 1, 2, 3, 
and 4 modes using a 1 Hz bandwidth with single point frequency samplings at each mode 


for a spatially complete beam. 
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VII. CONCLUSIONS / RECOMMENDATIONS 


A. SUMMARY 


Frequency Domain Structural Identification using the structural synthesis 
transformation was performed on a simulated experimental free-free beam model. The 
identification was performed for both the spatially complete and the spatially incomplete 


case. 


e SST based structural identification provides an exact solution for the 


identification of FE modeling errors given spatially complete data. 


e For spatially incomplete systems SST based structural identification provides 
a frequency dependent solution which is not suitable for finite element modeling error 


correcting. 


e Using both matrix and integral based techniques the SST can provide single 
mode solutions which are frequency independent correction matrices AK, AM, AC 
which approximately corrects the reduced FE model within a frequency bandwidth of 


the experimental system mode under consideration. 


B. CONCLUSIONS 


This thesis has clearly demonstrated that for spatially incomplete systems, single 
mode frequency independent solutions can be found which correct the reduced finite 
element model in a neighborhood of a given mode of the experimental system. It has also 
been shown that the concept of a single mode solution can be expanded to that of multiple 
mode solutions which are frequency independent correction matrices which approximately 
corrects the reduced FE model throughout a frequency bandwidth which includes more 


than 1 mode of the experimental system. 


C. RECOMMENDATION 


The purpose of this thesis was to find frequency independent multiple mode 


solutions which could be used to approximately correct reduced finite element models. 
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Although satisfactory results were obtained, investigation is still required in the following 


areas: 


eDetermine the relationship between errors in a full FE model and those of an 


associated reduced FE model. 


eDetermine a method to ‘pullback’ reduced FE model correction to full order 


FE model corrections. 


eInvestigate use of the integral formulation of reference (5) in analysis of the O- 
set system. 
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APPENDIX 


The following is a brief description of MATLAB routines employed in this thesis: 


*SST.M - Generates MAT files containing spatially complete or spatially incomplete 
experimental FRF, analytic FRF, localization matrix, and SST solutions. 


*PLOTSST.M - Uses files produced by SST.M to generate plots used in chapters 3 and 4. 


*CHAP3.M - Sets parameters and calls SST.M and PLOTSST.M to generate figures for 
chapter 3. 


*CHAP4.M - Sets parameters and calls SST.M and PLOTSST.M to generate figures for 
chapter 4. 


*CHAPS.M - Generate figures for chapter 5. 

*CHAP6.M - Generate figures for chapter 6. 

*SETUP.M - Generates MAT files for FE models. 

*BEAMMDL.M - Setup FE models. 

*FSTATIC.M - Perform static reduction of mass and stiffness matrices. 
*-FIRS TAM.M - Perform IRS reduction of mass and stiffness matrices. 
*FREQMODE.M - Returns FE model frequencies. 

*NDX3D.M - Indexes a series of 2-D matrices into a single 3-D matrix. 


*-INTSUB.M - Decomposes impedance matrix into mass, stiffness, and damping matrices using 
integral techniques. 


*MYTRAP.M - Computes integrals using trapezoidal method 
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SST.M 


MPWKLVHVHVHAMVVHVMHVLUwMHVUVUMUVHVMAVAMVVYNWMKVMUVUMNHV% 
Mw*“%%V%General setupw%%%%%%VHUNUUVKVMUVUVUwUMUMV%N% 
WWRwWLWHVV%VwMHVVAVUWUVWUMVVUVUWVUVVHLNYV“KV“VUWMVMNVV%N% 
clear,cle ;%%clear workspace. 

closeall ;Y%close any open figure windows 


jesqrt(-1) 
load sstconf ;%load A-set and O-set 
aset(oset)=[]; 
aset_size=length(aset); 
oset_size=length(oset); 
load beamdata 
if length(oset)== 0 
complete=1; 
else 
complete=0; 
end 


YYsetup plot labels~wW%%HUVMHVWHUVVHUVVWUHVLVMUMUVWVHUVVVVVAVVVVVVYVVVV%V%% 
if (oset==[])&(posofMasserr(1 }~=0)&(posofStifferr(1 }~=0) 
casename=(‘Spatially Complete Mass & Stiffness error’); 
elseif (oset==[])&(posofMasserr(1 }=0) 
casename=(‘Spatially Complete Mass error’); 
elseif (oset==[])&(posofStifferr(1)~=0) 
casename=(‘Spatially Complete Stiffness error); 
elseif (posofMasserr( 1 }>=0)&(posofStifferr(1)~=0) 
casename=(‘Spatially InComplete Mass & Stiffness error’), 
elseif (posofMasserr(1 }~=0) 
casename=(‘Spatially inComplete Mass error’); 
elseif (posofStifferr(1)}—=0) 
casename=(‘Spatially InComplete Stiffness error’); 
end 


%%eget reduced A set & O set frequenciessAKZ%UKKUVUWM%N%N%% 


kexto=k_anal(oset,oset); %% stifiness 
mexto=m_anal(oset,oset), %% and mass matrices 
if complete 
kstat=k_anal; 
mstat=m_anal; 
else 
if static == 0 
[kstat,:mstat]=fstatic(k_anal,m_anal,oset,aset);%% get reduced K & M 
elseif static == 
[kstat,mstat]=firs_tam(k_anal,m_anal,oset,aset);%% get reduced K & M 


else 
kstat=k_anal(aset,aset); 
mstat=m_anal(aset,aset); 
end 
end 


cstat=sqrt(-1)*struc_damping. *kstat; 
[u,lambdaa,c]=freqmode({k_anal,m_anal) ;%get freqs 
{u,lambdared,c]=freqmode(kstat,mstat) ;% 
[u,lambdax,c]=freqmode({k_exp,m_exp) 53% 
[u,lambdaexto,c]=freqmode(kexto,mexto),% 
omegaa=sqrt(lambdaa); 

omegax=sqrt(lambdax),; 

omegared=sqrt(lambdared); 
omegaexto=saqrt(lambdaexto),; 

fori=1:4 % this to delete fixed body modes 
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SST.M 


if real(omegax(1)) < 10“%(-3) 
omegax=oinegax(2:length(omegax)); 
end 
if real(omegaa(1)) < 10*(-3) 
omegaa=omegaa(2:length(omegaa)); 
end 
if complete==0 
if real(omegared(1)) < 10-3) 
omegared=omegared(2:length(omegared)); 
end 
if real(omegaexto(1)) < 10%(-3) 
omegaexlo=omegaexto(2:length(omegaexto)); 
end 
end 
end 


V*w“VY%%tind true mass and stiffness errors%w%%%%V%%KUMMUUMMNMUVUU%% 
true_stiffness=k_exp-k_anal:; 
stiftness_cset=find(diag(true_stiffness)); 
true_mass=m_exp-m_anal; 
mass_cset=find(diag(true_imass)); 
true_damping=c_exp-c_anal; . 
damping _cset=find(diag(true_damping)); 
if length(mass_cset) > 0 

for 1=1:length(mass_cset) 

x=find(stiffness_cset==mass_cset(1)); 


if ~length(x) > 0) 
stiffhess_cset=[stiffness_cset, mass_cset(1)]; 
end 
end 
end 


if length(damping_cset) > 0 
for 1=1:leugth(damping_cset) 
x=find(stiffness_cset==damping_cset(1)); 


if ~length(x) > 0) 
stiffness cset=[stiffness_cset; damping cset(i)]; 
end 
end 
end 


true_cset=sort(stiffness_cset) 

save true_errs true_cset true_mass true_stiffness true_damping 

clear true_mass true_stiffness true_damping 

clear area ee11 force g gc gcx gk gkx gm gmx goble goblcx goblk goblkx 

clear goblm gobImx ke lambdaa lambdaexto lambdared lambdax Iumpdamp 
clear lumpmass lumpspring me pho posofDamperr posofMasserr posof{Stflerr 
clear stifiness_cset 


%*%setup frequency range for analysis*H%~UVKKURAUVVKwVKHVHKHUVUVMVYUVUKHVYVAVL“LV“VVN% 
freqtop=omegax(highmode}+. I *omegax(highmode); 

freqbottom=omegax(lowmode)-. 1 *omegax(lowmode), 

freqbottom=10*2*pi; start at 10 hz 

freqtop=1200*2*pi; %end at 1200 hz 

%*Y%number of points to plot 


numpoints=400;%%%%%SET NUMBER OF POINTS TO USE FOR THE SIMULATION 
fineness=numpoints; 


w= freqbottom:(freqtop-freqbottom)(numpoints-1 ):freqtop; 


YWVLocate C setHwMwM%wHUW%%MM%%%M%MVMMMUMUMUMUVUVWU%MVURWUWMUMM%MVM%U%MUVMM%UUVUWVY 
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150 


155 


160 


165 


170 


175 


180 


keepgoing="'n' 
wireq=omegax(lowmode) ; 
nextmode=lownode; 

tol=1 i 
defaultstepsize=. 1; 

while (keepgoiug=='n’) 


z_anal_red=kstat+j*wfreq*cstat-wireq*2 *mstat ; 


h_anal_red=inv(z_anal_red) 


? 


z_exp=k_exptj* wfreq*c_ Srraoelee in EXD: 


h_exp=inv(z_exp) 
h_exp=h_exp(aset,aset) 
L=z_aual_red*(h_anal_red-h_exp 
ci=find(abs(diag(L))>tol) : 
temp_cset_size=length(c1); 
if length(c1) >= 1 
cset_size=length(c1) 
figure(1) 
plot(aset, abs(diag(L))) 
hold ou 


? 


yz. anal_red 


° 
> 


* 
? 


plot(aset(c1),zeros( 1 ,length(c1)),'x') 


tallest=max(abs(diag(L))) 
title((TOL = ‘num2str(tol 
ylabel(‘Ibf/in') 
if complete 

xlabel(‘DOF’) 
else 

xlabel(‘ASET DOF’) 
end 
hold off 


° 
? 
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),' Cset size = ‘int2str(cset_size),' Omega = ‘nun2str(wtreq/(2*p1)),' Hz'}) 


k=meuu( 'Clioose an action ‘,... 
‘Increase Tolerance ‘.... 
‘Decrease Tolerauce ‘,... 
‘Increase frequency '.... 
‘Decrease frequency ‘... 
' Previous mode .... 


' Next Mode 


i] 
peee 


‘Set Print switch on’... 
‘Increase stepsize ‘.... 
' Decrease stepsize ',... 


Go 
if k== 


?; 


stepsize=defaultstepsize; 


ctr=0; 


while length(ci) >= temp_cset_size 


ctr=ctr+1; 
if ctr==10 


stepsize=stepsize* 10; 


ctr=0 ; 
end 


tol=min(abs(tol+stepsize*(tol)),tallest) ; 


ci=find(abs(diag 
if lengthi(c1) >= 1 


(L)ptol) : 


cset_size=length(c1) : 


figure(1) 


plot(aset, abs(diag(L))) ; 


hold on 


plot(aset(c1),zeros(1,length(c1)),'x') 
tallest=max(abs(diag(L))) —; 


title({'TOL = 
ylabel(‘Ibf/in') 


‘num2str(tol),’ Cset size = 


"int2str(cset_size),’ Omega = 'num2str(wireq/(2*p1)),’ Hz']) 
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if complete 
xlabel("DOF’) 
else 
xlabel(‘ASET DOF') 
end 
hold off 
end 
end 
elseif k== 
stepsize=defaullstepsize; 
ctr=0; 
while length(ci) <= temp_cset_size 
ctr=ctr+l; 
if ctr==10; 
stepsize=stepsize*2; 
ctr=0; 
end 
tol=max(abs(tol-stepsize*(tol)),.00001) ; 
ci=find(abs(diag(L))>tol) : 
if length(ci) >= | 
cset_size=length(c1) 2 
figure(1) 
plot(aset, abs(diag(L))) : 
hold on 
plot(aset(ci),zeros(1,length(ci)),'x') 
tallest=max(abs(diag(L))) —; 
title([TOL = 
Hz']) 
ylabel(‘Tbf/in') 
if complete 
xlabel("‘DOF’) 
else 
xlabel(‘ASET DOF') 
end 
hold off 
end 
end 
elseif k==3 
while length(c1) == temp_cset_size 


"num2str(tol),’ Cset size = 


‘imt2str(cset_size),’ Omega = ‘nun2str(wtreq/(2*p1)), 


wireq=min(wfreqt.0000 1 *( freqtop-wireq),freqtop), 
z_anal_red=kstat+j*w(count)*cstat-wireq’2*imstat ; 


h_anal_red=inv(z_anal_red) : 


Z_exp=k_expty*wfreq*c_ Seed 2 Dene 


h_exp=inv(z_exp) 
h_exp=h_exp(aset,aset) 


L=z_anal_red*(h_anal_red-h bee anal_red 


tallest=max(abs(diag(L))) —; 
tol=nn(tol+.00 1 *tallest,tallest) ; 
ci=find(abs(diag eee : 

cset_size=leneth(ci) g 

figure(1) 

plot(aset, abs(diag(L))) : 

hold on 
plot(aset(c1),zeros(1,length(c1)),'x') 


title([{ TOL = 'num2str(tol),' Cset size = ',int2str(cset_size),’ Omega 


ylabel(‘Ibf/in') 
if complete 
xlabel(‘DOF') 
else 
xlabel((ASET DOF') 
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end 
hold off 
end 
elseif k==4 
wireq=max(wtreq-.05*(wfreq-freqbottom), freqbottom); 
elseif k==5 
nextmode=max(lowmode,nextinode-1 ); 
wtfreq=omegax(nextmode); 
elseif k==6 
nextmode=mun(highmode,nextmodet 1 ), 
wfreq=omegax(nextmode); 


elseif k== 
pswitch='y’ : 
elseif k==8 
defaultstepsize=defaultstepsize/10; 
elseif k== 
defaultstepsize=defaultstepsize* 10; 
else 
keepgoing='y’ ; 
cset=aset(ci); 
cset_rel=ci : 
end 
else 
tol=tol*.9; 
end 
end 
close(1) 


save extrtset cset cset_rel 
ARUV%%M%2%M%MMMUMMHAWMUMMNMMMUMMNMHNMAWANMMMUMUMUMUMUVMMMVMUVWMAHAUVUUVMUMUWUMUNMUMUVHUNMMV%V% 
YH*% 
YUVHAwVVMAwVMUVMM*wUWVMVWMUMVM““VFORCE THE CSETHW%%%%%M%%%%%M%M%MM%MNMM%MM%MMMM%MMMWHVVO% 
Ycset=[7 9 11 13] %%%%%%%%THE INCOMPLETE CASE CSET%%*%%%%%%%M%N%%% 
%cset_rel=[4 5 6 7]%%%%%%%WM%WMVWHWUUMVVM*WMUMYWMUHUVVUWHUWVWVHVWMUMUM%%MMM%M%%MM%M%%P*7MM%N% 
YMRHRH%%%%%%%W%W%%MMUW%W%W%%MM%%MM%MM%%MU%MU%-%MMVWMUVYWVVW%M%MWwVYV*w%MUMVMUVUVVWVW%VVWVVWWV%N 
YVV% 
%cset=[4 56789 10 11J%WW%%%%%R%W%M%%H%MYVTHE COMPLETE CASE CSET%%% 
Ycset_rel=[(4 56789 10 LLJ%R*Z%X*XwRKKHRwWKVUwVUWRwVUwHUUHUMRWVWHUU%MMUMMV%%%V%VV%%V% 
YMH%UH%%%%%%%%%%MU%UHU%MU%MM%MHMM%%%MMMMM%%%NMMUHMMMHMMHMMVUMMV%%MMMVM% 
MRwH% 
cset_size=length(cset); 
save LL 
clear L 
clear kexta kexto mexta mexto 
pack 
skyfull=zeros(numdof,numdof); 
count=1 ; 
for index] =1:numdof 
for index2=index 1 :numdof 
skyfull(index]1 ,index2)}=count; 
skyfull(index2,index 1 }=count; 


count=countt 1; 
end 
end 
skyred=zeros(aset_size,aset_size); 
count=]  ; 


for index1=l:aset_size 
for index2=index] :aset_size 
skyred(index] ,index2)=count; 
skyred(midex2,index 1 )=count, 
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count=count+ 1; 


end 
e1id 
skycset=zeros(cset_size,cset_size); 
count=1 ; 


for index1=1 :cset_size 
for index2=index 1 :cset_size 
skycset(index 1 ,index2)=count; 
skycset(index2,index 1 )=count; 
count=count+ 1; 
e1d 
end 


full_holder=zeros(1,(numdof*(numdof+1 ))/2); 
red_holder=zeros(1,(aset_size*(aset_size+1))/2); 
cset_holder=zeros(1,(cset_size*(cset_size+| ))/2), 


if meters 
waitbar_handle=waitbar(0,'Computing Experimental FRF'), 
else 
disp(‘Getting experimental FRF") 
end 
H_EXP={]; 
for count=1:numpoints 
if meters 
waitbar(count/numpoints); 
end 
Z_exp=k_exptj*w(count)*c_exp-w(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex= 1; 
for index1=1:aset_size 
for index2=1ndex 1 :aset_size 
red_holder(skyindex)=h_exp(index]1 ,index2), 
skyindex=skyindex+ 1; 
end 
end 
H_EXP=[H_EXP red_holder’}; 
end 
if meters 
close(waitbar_handle) 
end 
save H_EXP H_EXP 
clear H_EXP 
pack 


if meters 
waitbar_handle=waitbar(0,'Computing Expernnental lLmpedence’), 
else 
disp(‘Getting experimental Impedence') 
end 
Z_EXP=[]; 
for count=1:numpoints 
if meters 
waitbar(count/numpoints), 
end 
z_exp=k_exptj*w(count)*c_exp-w(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
Z_exp=mv(h_exp); 
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skyindex=1 ; 
for index]=] :aset_size 
for index2=index ]:aset_size 
ted_holder(skyindex)=z_exp(index 1 ,index2); 
skyindex=skyindex+1; 
end 
end 
Z_EXP=[Z_EXP red_holder); 
end 
if meters 
close(waitbar_handle) 
e1d 
save ZA EX? Z EXP 
clear Za XP 
pack 


if meters 
waitbar_handle=waitbar(0,'Computing Reduced Analytic FRF'), 
else 
disp(‘Getting Reduced FRF’) 
end 
H_ANAL RED=[]; 
for count=1] :umpoints 
if meters 
waitbar(counUnumpoints),; 
end 
z_aual_red=kstat+j*w(count)*cstat-w(count)*2 *imstat; 
h_anal_red=inv(z_anal_red), 
skyindex=1 ; 
for uidex1=1:aset_size 
for index2=index] :aset_size 
red_holder(skyindex)=h_anal_red(index 1 ,index2), 
skyindex=skyuidex+ 1; 
end 
end 
H_ANAL_RED=[H_ANAL_RED red_holder’]; 
end 
if meters 
close(waitbar_handle) 
end 
save H ANAL RH_ANAL_ RED 
clear Z_ ANAL RED H_ANAL_ RED kexta kexto mexta mexto 
pack 


%% coinpute and save L matnx diagonal as a function of frequency%%%w%%%%%% 
if ineters 
waitbar_handle=waitbar(0,'Computing L Diagonals’), 
else 
disp(‘Getting L Matrix’) 
end 
L_DIAGS=[]; 
for count= | :numpoints 
if ineters 
waitbar(count/numpoints ); 
end 
z_anal=k_anal+j*w(count)*c_anal-w(count)*2*m_anal; 
z_anal_red=Kstat+j*w(count)*cstat-w(count)*2 *mstat; 
h_ana]_red=inv(z_anal_red); 
Z_exp=k_exptj*c_exp-w(count)*2*m_exp; 
h_exp=inv(z_exp),; 
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h_exp=h_exp(aset,aset); 
L=z_anal_red*(h_anal_red-h_exp)*z_anal_red; 
temp_1_diags=diag(L); 
L_DIAGS=(L_DIAGS temp_1_diags(:)]; 

eud 

if meters 
close(waitbar_handle) 

end 

save L DIAGS L_DIAGS 

clear L_ DIAGS z_anal temp _1 diags L 

clear z_exp z_anal_reduv h_exp h_anal_red 


pack 


if meters 
waitbar_handle=waitbar(0,'Computing Z Analytic’); 
else ~ 
disp(‘Getting Z anal’) 
end 
Z_ANAL=[]; 
for count=1:numpoints 
if meters 
waitbar(count/numpoints); 
end 
z_anal=k_anal+j*w(count)*c_anal-w(count)*2*m_anal; 
skyindex=1_ ; 
for index 1=1:numdof 
for index2=index | :numdof 
full_holder(skyindex)=z_anal(index1,imdex2); 
skyindex=skyindex+ 1; 
end 
end 
Z_ ANAL=[Z_ANAL full_holder'’; 
end 
if meters 
close(waitbar_handle) 
end 
save Z_ANAL Z_ANAL 
clear Z ANAL z_anal 
pack 
if meters 
waitbar_handle=waitbar(0,'Computing Reduced Z Analytic’); 
else 
disp(‘Getting Reduced Z anal’) 
end 
Z_ANAL RED={]J; 
for count=]:numpoints 
if meters 
waitbar(counU/numpoints); 
end ; 
z_anal_red=kstat+j* w(count)*cstat-\w(count)*2* mstat; 
skyindex=1_ ; 
for index1=1:aset_size 
for index2=index! :aset_size 
red_holder(skyindex)=z_anal_red(index 1 ,index2); 
skyindex=skyindex+ 1; 
end 
end 
Z_ ANAL _RED=[Z ANAL RED red_holder]; 
end 
if meters 
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close(waitbar_handle) 
end 
save Z ANAL _RZ_ANAL_ RED 
clear Z ANAL RED z_anal_red 
485 pack 


YVM*VYcompute DZ%%%V%%%VYVUVWUVVUUVUMUVUVUVUVUV%VVUVVVVVVVUVVVUVUVVVUVUVVVUVHUYUVVV% 
if meters 
waitbar_handle=waitbar(0,'Computing DZ...'); 
490 else 
disp(‘Getting DZ’) 
end 
DzZ=(); 
for count=1:numpomts 
495 if meters 
waitbar(count/numpoints); 
end 
z_anal_red=kstatt+j*w(count)*cstat-w(count)*2 *mstat; 
h_anal_red=mnv(z_anal_red); 
500 Z_exp=k_exptj*w(count)*c_exp-w(count)*2*m_ exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
505 =inv(inv(inv(hacc )*(hacc-hxec)* inv(hacc))-hacc); 
skyindex=1_ ; 
for index1=1 :cset_size 
for index2=index | :cset_size 
cset_holder(skyindex)=dz(index] ,index2),; 


510 skyindex=skyindex+] ; 
end 
end 
DZ=[DZ cset_holder’}; 
end 


SD if meters 
close(waitbar_handle) 
end 
save DZ DZ 
clear hace hxee dz h_exp z_anal_red h_anal_red z_exp 
920 pack 


%*~%Ydecompose DZ into DM, DK, & DC%V%%W%w%%UR%wMNMMMUM%UMUMUMUVUUMVUVVMAVVMVVUVVAUWAAVV%V% 
if meters 
waitbar_handle=waitbar(0,'Decomposing DZ...'); 
525 else 
disp(‘Decomposing DZ’) 
end 
DK=(]; 
DM=(); 
530 DC={]; 
dz1=zeros(cset_size,cset_size); 
=zeros(Cset_size,cset_ size); 
dz3=zeros(cset_size,cset_size); 
for i=] :numpoints-2 
535 if meters 
waitbar(i/(numpoints-2 )) 
end 
dztemp=DZ(ndx3d([1 cset_size*(cset_size+] )/2 numpoints],1,]:cset_size*(cset_size+1)/2,1)), 
skyindex=1, 
540 for udex1=1:cset_size 
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for index2=index 1 :cset_size 
dz1(index 1 ,index2)=dztemp(skyindex ); 
dz1(index2,index 1 }=dztemp(skyindex), 
skyindex=skyindex+1; 
end 
end 
dztemp=DZ(ndx3d({1 cset_size*(cset_size+1)/2 nwnpoints], 1, 1:cset_size*(cset_size+1)/2,i1+1)), 
skyindex=1; 
for index1=1:cset_size 
for index2=index]:cset_size 
dz2(index 1 ,index2)=dztemp(skyindex), 
dz2(index?2, index 1 }=dztemp(skyindex),; 
skyindex=skyindex+ 1; 
end 
end 
dztemp=DZ(ndx3d([1 cset_size*(cset_size+1)/2 numpoints], 1, 1:cset_size*(cset_size+1)/2,i+2)), 
skyindex=1; 
for index1=1:cset_size 
for index2=index I :cset_size 
dz3(index 1 ,index2)=dztemp(skyindex); 
dz3(index2,index 1 )=dztemp(skyindex); 
skyindex=skyindex+ 1; 
end 
end 
temp=[eye(cset_size) -w(i)*2*eye(cset_size) +j*w(i)*eye(cset_size); 
eye(cset_size) -w(it+]1 )*2*eye(cset_size) +)*w(it+1 )*eye(cset_size); 
eye(cset_size) -w(1+2)*2*eye(cset_size) +j*w(it2)*eye(cset_size)]\[dz1;dz2;dz3], 
ktemp=temp(1:cset_size,:); 
mtemp=temp(cset_size+]:2*cset_size,:); 
ctemp=temp(2*cset_size+1:3*cset_size,:), 


skyindex=1 —; 
for index1=1:cset_size 
for index2=index] :cset_size 
cset_holder(skyindex)=ktemp(index] ,index2), 
skyindex=skyindex+1; 
end 
end 
DK=([DK cset_holder(:)]; 


skyindex=1 —; 
for index]=1:cset_size 
for mdex2=index]:cset_size 
cset_holder(skyindex)=mtemp(index 1 ,index2), 
skyindex=skyindex+1; 
end 
end 
DM=[DM cset_holder(:)}: 


skyindex=] 5; 
for index1=1:cset_size 
for index2=index] :cset_size 
cset_holder(skyindex)=ctemp(index 1 ,index2), 
skyindex=skyindex+ 1; 
end 
end 
DC=[DC cset_holder(:)}, 


end 
if meters 


close{waitbar_handle) 
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end 
save DM DM 
save DK DK 
save DC DC 

605 clear DC DK DM 
save exp c_exp k_exp m_exp 
clear c_exp k_exp m_exp 
save stat kstat mstat cstat 
clear kstat mstat cstat 

610 save anal k_aual m_anal c_anal 
clear k_anal m_anal c_anal 
pack 


if complete 
GTS load Z_ANAL 
Z=Z_ ANAL; 
clear Z_ ANAL 
else 
load Z_ ANAL R 
620 Z=Z_ANAL_RED; 
clear Z ANAL RED 
end 
CORRHA={]; 
tempza=zeros(aset_size,aset_ size), 
625 tempz=zeros(cset_size,cset_size), 
if meters 
waitbar_handle=waitbar(0, Installing DZ...'), 
else 
disp(‘Installing DZ’) 
630 end 
for i= 1:numpoints 
if meters 
waitbar(i/numpoints) 
end 
635 ztemp=Z(ndx3d([1 aset_size*(aset_size+1)/2 (numpouits-(1-1 ))],1,]:aset_size*(aset_size+1)/2,1)), 
skyindex=1; 
for index 1=]:aset_size 
for index2=index|]:aset_size 
tempza(index 1 ,index2)=ztemp(skyindex); 


640 tempza(index2,index 1 )=ztemp(skyindex); 
skyindex=skyindex+1, 
end 
end 
Z:1 (I; 
645 ztemp=DZ(ndx3d([1 cset_size*(cset_size+])/2 (numpoints-(1-1))],1, 1:cset_size*(cset_size+1)/2,1)); 


skyindex= 1, 
for index 1=1:cset_size 
for index2=1ndex]:cset_size 
tempz(index | ,index2 )=ztemp(skyindex); 
650 tempz(index2,index 1 }=ztemp(skyindex), 
skyindex=skyindex+ 1; 
end 
end 
tempza(cset_rel,cset_rel)=tempza(cset_rel,cset_rel}+tempz; 
Boo DZ¢:, 1 ={]; 
tempha=inv(tempza), 
CORRHA=[CORRHA tempha(:)}; 
end 
if ineters 
660 close(waitbar_ handle) 
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end 
clear DZ Z 
save CORRHA CORRHA 
clear CORRHA 
665 
load DK; 
load DM; 
load DC; 


670 load stat 
CORRHAD=[]; 
tempza=zeros(aset_size,aset_size); 
tempz=zeros(cset_size,cset_size); 
tempk=zeros(cset_size,cset_size); 
675 tempm=zeros(cset_size,cset_size); 
tempc=zeros(cset_size,cset_ size); 
if meters 
waitbar_handle=waitbar(0,‘Installing DK, DM, & DC...'); 
else 
680 disp(‘Installing DK, DM, DC’) 
end 
for 1=1:nwmnpoints-2 
if meters 
waitbar(i/(numpoints-2 )) 
685 end 
kcorrected=kstat; 
mcorrected=mstat; 
ccorrected=cstat; 


690 temp=DK(ndx3d({1 cset_size*(cset_size+1)/2 (numpoints-(2))],1,1:cset_size*(cset_size+1)/2,1)); 
skyindex= 1; 
for index1=1:cset_size 
for mdex2=1ndex1 :cset_size 
tempk(index 1 ,index2)=temp(skyindex); 
695 tempk(index2, index 1 )=temp(skyindex); 
skyindex=skymdex+1; 
end 
end 


700 DKc:,1}=[]; 
temp=DM(ndx3d([1 cset_size*(cset_size+1)/2 (numpouits-(2))],1,1:cset_size*(cset_size+1)/2,1)); 
skyuidex= 1; 
for index1=1:cset_size 
for index2=index 1 :cset_size 
705 tempin(index 1 ,mdex2 =temp(skyindex); 
tempm(index2,index 1 }=temp(skyindex); 
skyiidex=skyindex+1, 


end 
end 
710 
DM(:,1)=[]; 
temp=DC(ndx3d([1 cset_size*(cset_size+1)/2 (numpoints-(2))], I, l:cset_size*(cset_size+1)/2,1)); 
skyindex=1, 
for index1=1:cset_size 
715 for undex2=index] :cset_size 
tempc(index 1 ,index2 )=temp(skyindex); 
tenipe{index2, index] )=temp(skyindex), 
skymdex=skyindex+1, 
end 
420 end 
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DC(:,D=U; 
kcorrected(cset_rel,cset_rel)}=kcorrected(cset_rel,cset_rel}+tempk; 
mcorrected(cset_rel,cset_rel)=incorrected(cset_rel,cset_rel}+teinpin; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel}+teinpe; 


tempza=kcorrected+j* w(it 1 )*ccorrected-w(it+1)*2*mcorrected; 
tempha=inv(tempza),; 
CORRHAD=[CORRHAD templia(:)]; 
end 
if meters 
close(waitbar_handle) 
end 
clear tempk tempm tempc tempza tempha 
clear kcorrected mcorrected ccorrected 
save CORRHAD CORRHAD 
clear CORRHAD DK DM DC 


M%Vsave the 
Workspace%%%%%%%%%%%%M%MMMM%M%MMMMMMMMMUMMUVAHANMNVVMMMMVMUMMVMMM%%N% 
save INT 
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clear 
closeall 
load INT 
titles=0; 
if pswitch== 
whitebg(‘white') 
close 
end 
h=0; 
if complete 
fignum=1; 
else 
fignum=2 1; 
end 
%% write frequencies and system diag to diary file prt.outw»/w%nwnynwnwAU%wVRw%w%U%%MNMNM%%N% 
if exist(fig000. out’) 
delete fig000.out 
end 
diary fig000.out 
fprintf" — \n') 
fprint{(' \n') 
tprintf(‘A set\n') 
for index=1:length(aset) 
fprintf{({blanks(4-length(int2str(aset(index)))),int2str(aset(index))}) 
if rem(index, 12)==0 
fprintf(‘\n') 
end 
end 
fprintf{(\n‘) 
fprintf(‘\n') 
fprint{(‘O set\n') 
for index=1:length(oset) 
tpnntf{([blanks(4-length(int2str(oset(index)))),int2str(oset(index))]) 
if rem(index,12)==0 
fprintf{(‘\n') 
end 
end 
fprintf(‘\n') 
{printf(‘\n') 
fprintf((Computed C set\n') 
for index=1:length(cset) 
fprintf([blanks(4-length(int2str(cset(index)))),int2str(cset(index))]}) 
if rem(index,12)== 
fprintf(‘\n') 
end 
end 
fprintf(‘\n') 
tprintf(\n') 
fprint{(‘True C set\n’) 
for index=1:length(true_cset) 
fprintf([(blanks(4-length(int2str(true_cset(index)))),int2str(true_cset(index))]) 
if rem(index, 12)==0 
fprintf(‘\n') 
end 
end 
fprintf(\u') 
fprintf(\n') 
bigger=max(length(omegaa),length(omegax)); 
if complete 
allfreqs={ ... 
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[omegaa/(2*pi); zeros( bigger-length(omegaa), 1 )] ... 
{omegax/(2*pi), zeros(bigger-length(omegax), 1)] ... 

else 
allfreqs=[ ... 
{omegaa/(2* pi); zeros(bigger-length(omegaa),1) J... 
{omegax/(2*pi); zeros(bigger-length(omegax),1) _—j... 
[omegared/(2*pi); zeros(bigger-length(omegared),1) ] ... 
{omegaexto/(2*pi);, zeros(bigger-length(omegaexto), 1 )]... 


iF 
end 
{aaa bbb]=size(allfreqs); 
fprintf(‘System Frequencies (Hz)\n’) 
if complete 
fprintf( Anal Exp\n') 
bbb=2; 
else 
fpnntf’? Anal Exp Red Oset\n') 
bbb=4, 
end 
for index]=]:aaa 
prtline=[]; 
for index2=1:bbb 
if allfreqs(index] ,index2)== 
prifreq=" 
else 
prtfreq=[blanks(8-length(sprintf('%o4 g' allfreqs(index 1 ,index2)))).... 
sprintf('%o4g',allfreqs(index1,index2))}; 
end 
prtline=[prtline, prtfreq]; 
end 
prtline=[prtline, ‘\n’]; 
fprintf(prtline) 
end 
diary oft 
%dos(‘type figd00.out > Iptl &’); 
xx=omegax(find(omegaa(lowmode) <= omegax & omegax <= omegaa(highmode)))(2*p1),; 
yx 1=ones(length(xx), 1); 
aa=omegaa(find(omegaa(lowmode) <= omegaa & omegaa <= omegaa(highimode)))/(2 *pi); 
yal=ones(length(aa),1); 
red=omegared(find(omegaa(lowmode) <= omegared & omegared <= omegaa(highmode)))/(2 *p1); 
yred1=ones(length(red), 1); 
exto=omegaexto(find(omegaa(lowmode) <= omegaexto & omegaexto <= omegaa(highmode)))/(2* pi); 
yextol=ones(length(exto), 1); 
if complete 
x0=max([xx(length(xx)) aa(length(aa))]); 
else 
x0=max([xx(length(xx)) aa(length(aa)) red(length(red)) exto(length(exto))}); 
end 
clear omegax omegaa omegared omegaexto 
clear LZ ANAL c_anal c_exp conn h_exp h_anal_red k_anal k_exp 
clear kexta kstat kexto m_anal m_exp mexto mstat temp_l_diags 
if complete 
mid_index=round(length(cset)/2); 
else 
mid_index=round(length(cset)/2); 
end 
load H_ANAL_R 
load H_EXP 
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h=figure(h+1); 
ficnum=fignuint+| ; 
plot(w/(2*pi),log10(abs(H_ANAL_ RED(ndx3d([1 aset_size*(aset_size+])/2 numpoints].... 
1 ,skyred(cset_rel(imid_index),cset_rel(mid_index)),':')))),'r—‘... 
125 W/(2*pi),logiO(abs(H_EXP(ndx3d([1 aset_size*(aset_size+1)/2 nwmpoints],... 
1,skyred(cset_rel(inid_index),cset_rel(imid_index)),""')))),’b-.') 
V=anIS; 
ylabel(("H(‘,int2str(cset(nud_index)),',",mt2str(cset(mid_index)),") mw/lbt (LogiO of)'}) 
xlabel((Omega (Hz)') 
130 if titles 
title(‘Analytic FRF vs Experimental FRF') 
end 
hold on 
if abs(v(4}-v(3)) <= 1 
13S v2(1)=v(1); 
v2(2)=v(2); 
v2(3)=v(3}-. 1 *abs( v(4)); 
v2(4)=v(4)}+. 1 *abs(v(4)); 
axi1s(v2); 
140 end 
grid on 
hh=slegend(ANALYTIC', EXPERIMENTAL’); 
axes(hh); 
hold off 
145 if pswitch=='y' 
prtfig(fignum) 
delete(h) 
end 


150 clear H-ANAL RED 
clear H_EXP 


if complete 
h=figure(h+1); 
(bere, fignum=fignumt |; 
else 
h=figure(h+2); 
fignum=fignum+2; 
end 
160 load L 
plot(aset,abs(diag(L)),aset,abs(diag(L)),'*') 
ylabel(' L(DOF) Ibffin') 
if complete 
xlabel(‘DOF') 
165 else 
xlabel(‘ASET DOF') 
end 
erid on 
if titles | 
170 title(['Localization Matrix Diagonal at Omega = ‘num2str(wtreq/(2*p1)), ' Hz']); 
end 
if pswitch=='y' 
prtfig(fignum) 
delete(h) 
ro end 
clear L 


h=ligure(h+1); 


fignum=fignum+ 1; 
180 load L_DIAGS 
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mesh(w/(2*pi),aset,log10(abs(L_DIAGS))) 
if titles 

title(‘Freque1icy Dependence of Localization Matrix Diagouals') 
end 
xlabel(‘Omega (Hz)') 
ylabel(‘DOF’) 
zlabel('L(DOF, omega) Ibf/in (Log10 of)') 
grid on 
if pswitch== 

prtfig(fignum) 

delete(h) 
end 


h=figure(h+1 ); 
fionum=fignum+t1; 
subplot(2,1,1) 
plot(w/(2*p1),log 10(L_DIAGS(ndx3d([aset_size 1 numpoints],cset_rel(mid_index+2),1,"')))) 
ylabel([Error Coord L(’,int2str(cset(imid_index+2)),',"int2str(cset(mid_index+2)),")']) 
if titles 
title(‘Ibf/in (Log10 of) ); 
end 
grid on 
vl=axis; 
subplot(2,1,2) 


if length(cset_rel) < length(aset) 
holdset=1:length(aset) ; 
holdset(cset_rel)=[] ; 
index=round(length(holdset)/2), 


plot(w/(2*pi),log 1O(L_DIAGS(ndx3d({aset_size 1 numpoints],holdset(index), 1,':')))) 
ylabel(['Non-Error Coord L(’,int2str(aset(holdset(index))),',",int2str(aset(holdset(index))),')]) 
if titles 
title(Frequency Dependence of L Matnx Non-Error set Diagonal Elements' ) 
end 
v=axis; 
hold on 
v2(1}=v(1); 
v2(2)}=W(2); 
v2(3}=v(3}+((v(4 )-v(3))*.5 Hv1(4)-v1(3))*.5 ; 
v2(4)=v(3 H{(v(4)-v(3))*.5 #(v1(4)-v1(3))*.5; 
axis(v2); 
grid on 
end 
xlabel(‘(Omega (Hz)) 
hold off 
if pswitch=='y' 
prtfig(fignum) 
delete(h) 
end 
clear L_DIAGS 


load Z_ANAL_R 

load Z_EXP 

h=figure(h+ 1); 

fignumn=fignum+ 1; 

plot(w/(2*pi),log10(Z_ANAL_RED(ndx3d({1 aset_size*(aset_size+1)/2 

muupoints], 1 ,skyred(cset_rel(mid_index),cset_rel(inid_idex)),""))),'r—-',... 
W(2*pi),logl0(Z_EXP(ndx3d({1 aset_size*(aset_size+1)/2 

numpoints], | skyred(cset_rel(mid_index),cset_rel(mid_index)),"'))),'b-.’) 
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xlabel(‘(Omega’*2 (Hz)') 
ylabel([(Z(’,int2str(cset(mid_index)),",’int2str(cset(mid_iudex)),") Ibf/in (Log10 of)'}) 
if titles 
title(‘Analytic Impedance vs Experunental Impedance’) 
end 
Vv=anis; 
hold on 
if abs(v(4)-v(3)) <= 1 
v2(1)=w(1); 
v2(2)=v(2); 
v2(3)=v(3)-.1 *abs(v(4)); 
v2(4)}=v(4)+. 1 *abs(v(4)); 
axis(v2); 
end 
erid on 
if complete== 0 
hh=slegend(‘ANALYTIC IMPEDANCE','EXPERIMENTAL IMPEDANCE’), 
else 
hh=slegend(ANALYTIC IMPEDANCE'’,EXPERIMENTAL IMPEDANCE), 
end 
axes(hh) : 
hold off 
if pswitch=='y' 
prtfig(fignum) 
delete(h) 
end 
clear Z_ ANAL RED Z EXP 


load true_errs 
load DK 
h=figure(h+1) ; 
fignum=fignum+] ; 
if complete 
plot(w(1:numpoints-2 )/(2 *pi), DK(ndx3d([{1 cset_size*(cset_size+1)/2 numpoints-2],1,skycset(mid_index,mid_index),'")),1T— 
W(1:numpoints-2)/(2*pi),true_stiffness(cset(imid_index),cset(mid_index)).*ones(1,numpoints-2),'b-.') 
ylabel([(DK(‘,int2str(cset(mid_index)),’,",int2str(cset(mid_index)),") 1bf/in']) 
xlabel(‘Omega (Hz)’) 
else 
plot(w( 1 :numpoints-2 (2 *pi),log 10(abs(DK(ndx3d({1 cset_size*(cset_size+1] /2 numpoints- 
2],1,skycset(mid_index,mid_index),'"')))),T—~'... 
W(I:numpoints-2)/(2*p1),log10((true_stiff{ness(cset(mid_index),cset(inid_index)) == 0)+... 
true_stiffness(cset(muid_index),cset(mid_index))).*ones(1,numpouits-2),’b-.") 
ylabel([‘DK(‘,int2str(cset(imid_index)),',’,int2str(cset(mid_index)),") 1bf/in (Log10 of)']) 
xlabel(‘Omega (Hz)') 
end 
V=axis 5 
if titles 
title(‘Computed Stiffness Error vs True Stiffness Error’) 
end 
v2(1)=v(1); 
v2(2)=v(2); 
if complete 
if abs(v(4)-v(3)) < 1 
v2(3)=v(3)-100*abs(v(4)-v(3)); 


v2(4}=v(4)+100* abs(v(4)-v(3)); 
v(3)=v2(3); 
v(4)=v2(4); 

end 
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if abs(true_stiffness(cset(mid_index),cset(mid_index)) - v(3)) < .25*abs(v(4)-v(3)) 
v2(3)=v(3)-.5*abs(v(4}-v(3)) ; 
v2(4)=v(4); 
elseif abs(true_stiffness(cset(mid_index),cset(mid_index)) - v(4)) < .25*abs(v(4)-v(3)) 
v2(4)=v(4 }+.5 *abs(v(4)-v(3)) ; 
v2(3)=v(3); 
else 
v2(3)=v(3); 
v2(4)=v(4), 
end 
else 
if abs(log 1 0((true_stiffness(cset(mid_index),cset(mid_index)) == 0)-... 
true_stiffness(cset(mid_index),cset(mid_index))) - v(3)) < .25*abs(v(4)}-v(3)) 
v2(3)=v(3)-.5*abs(v(4)-v(3)) ; 
v2(4)-v(4); 
elseif abs(log10((true_stiffness(cset(mid_index),cset(mid_index)) == 0}... 
true_stiffiress(cset(mid_index),cset(mid_index))) - v(4)) < .25*abs(v(4)-v(3)) 
v2(4)=v(4}+.5 *abs(v(4}-v(3)) ; 
v2(3)=v(3); 
else 
v2(3)=v(3), 
v2(4)=v(4); 
end 
end 
aX1S(v2) 
end on 
hh=slegend((COMPUTED STIFFNESS ERROR',TRUE STIFFNESS ERROR’), 
axes(hh) z 
if pswitch=='y' 
prtfig(fignum) 
delete(h) 
end 
clear DK 
load DM 
h=tigure(h+1) ; 
fignum=fignumt] ; 
if complete 
plot(w(1 :numpoints-2)/(2*p1), DM(ndx3d({1 cset_size*(cset_size+1)/2 numpoints-2],1,skycset(mid_index,mid_index),"")),'’r— 


4 
gree 


w(1:numpoints-2)/(2*pi),true_imass(cset(nud_index),cset(mid_index)).*ones(1,numpoints-2),'b-.') 
ylabel({DM(‘,int2str(cset(mid_index)),',int2str(cset(mid_index)),') lbf-sec*2/in'}) 
xlabel(‘Omega (Hz)') 
else 
plot(w( 1 :numpoints-2)/(2* pi), log 10(abs(DM(ndx3d([1 cset_size*(cset_sizet+1)/2 numpomnts- 
2},1,skycset(mid_mdex,nud_index),""')))),'r—".... 
W(1:numpoints-2 /(2 *pr),log 10(... 
& 
true_mass(cset(mid_index),cset(mid_index)) == 
}ttrue_mass(cset(mid_index),cset(mid_index))... 
).*ones(1,numpoints-2),’b-.') 
ylabel({'DM(,int2str(cset(mid_index)),',',int2str(cset(mid_index)),') lbf-sec*2/in (Log10 of)’)) 
xlabel(‘Omega (Hz)') 
end 
Vv=axis; 
if titles 
title((Computed Mass Error vs True Mass Error’) 
end 
hold on 
v2(1)=v(1); 
v2(2)=v(2),; 
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if complete 

if abs(v(4)-v(3)) < 1 
v2(3)=v(3)}-10000*abs(v(4)-v(3)); 
v2(4}=w(4)+10000* abs(v(4)-v(3)); 
v(3)=v2(3), 
v(4)=v2(4), 

end 

if abs(true_mass(cset(mid_index),cset(mid_tndex)) - v(3)) < .25*abs(v(4)-v(3)) 
v2(3)=v(3}-.5*abs(v(4)-v(3)); 
v2(4)=W(4): 

elseif abs(true_mass(cset(nud_index),cset(mid_index)) - v(4)) < .25*abs(v(4)}-v(3)) 
v2(4)=Vv(4)+.5*abs(v(4)-v(3)); 
v2(3)=v(3); 

else 
v2(3)=v(3); 
v2(4)=w(4): 

end 

else 

if abs(log10((true_mass(cset(mid_index),cset(muid_index)) == 0}... 
true_mass(cset(mid_index),cset(mid_index))) - v(3)) < .25*abs(v(4}-v(3)) 
v2(3)=v(3}-.5*abs(v(4)-v(3)) ; 
v2(4)=v(4); 


elseif abs(log10((true_mass(cset(nud_index),cset(mid_index)) == 0)-... 


true_mass(cset(imid_index),cset(mid_index))) - v(4)) < .25*abs(v(4)-v(3)) 
v2(4)=v(4)+.5 *abs(v(4)-v(3)) ; 
v2(3}-v(3); 
else 
v2(4)=v(4), 
v2(3)=v(3); 
end 
end 
axis(v2) 
grid on 
hh=slegend’ COMPUTED MASS ERROR',TRUE MASS ERROR’), 
axes(hh) ; 
if pswitch=="y' 
prtfig(fignum ) 
delete{h) 
end 
clear DM 


load DC 
h=figure(h+1 ), 
fignum=fignum+ 1; 
if complete 
plot(w(1 :numpoints-2)/(2*pi),imag(DC(ndx3d([1 cset_size*(cset_sizet+1)/2 numpounts- 
2], 1,skycset(inid_index,mid_index),'"))),'r—"'... 
w(1:numpomits-2)/(2*p1),imag(true_damping(cset(imid_index),cset(mid_index)).*ones(1,numpoints-2)),'b-.') 
ylabel({"DC(’,int2str(cset(mid_index)),’,’,int2str(cset(nmud_index)),") lbf-sec/in’}) 
xlabel((Omega (Hz)') 
else 
plot(w(1 snumpoints-2)/(2*pi),log10(abs( DC(ndx3d({1 cset_size*(cset_sizet+1)/2 numpoiits- 
2],1,skycset(imid_index,mid_index),"’)))),'r—... 
W(1 :numpoints-2)/(2*pi),log10(abs¢... 
fe 
true_damping(cset(mid_index),cset(mid_index)) == 0 ... 
}ttrue_damping(cset(mid_index),cset(mid_index))... 
).*ones(1 ,numpoints-2)),"b-.’) 
ylabel(['DC(‘,int2str(cset(mid_index)),’,’,int2str(cset(mid_index)),') Ibf-sec/in (Log10 of)'}) 
xlabel(Omega (Hz)’) 
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end 
v=anxis; 
if titles 
title(‘Computed Damping Error vs True Damping Error’) 
end 
v2(1)=w(1); 
v2(2)=v(2); 
if complete 
if abs(v(4)-v(3)) < 1 
v2(3)=v(3)-100 *abs( v(4)-v(3)); 
v2(4)}=v(4)}+100*abs(v(4)-v(3)); 
v(3)=v2(3); 
v(4)=v2(4); 
end 
if abs(true_damping(cset(mid_iidex),cset(mid_index)) - v(3)) < .25*abs(v(4)-v(3 )) 
v2(3)=v(3)-.25 *abs( v(4}-v(3)) ; 
v2(4)=v(4); 
elseif abs(true_damping(cset(mid_index),cset(imid_index)) - v(4)) < .25*abs(v(4)-vw(3)) 
v2(4)=v(4}+.25*abs(v(4)-v(3)) ; 
v2(3)}=v(3); 
else 
v2(3)}=v(3); 
v2(4)=v(4); 
end 
else 
if abs(log1 0((true_damping(cset(mid_index),cset(mid_index)) == 0}... 
true_damping(cset(mid_index),cset(imid_index))) - v(3)) < .1 *abs(v(4)-v(3)) 
v2(3)=(v(3)-.25 *abs(v(4 )-v(3))) 5 
v2(4)=W(4), 
elseif abs(log10((true_damping(cset(mid_index),cset(nud_index)) == O}... 
true_damping(cset(mid_index),cset(mid_index))) - v(4)) < .1 *abs(v(4)-v(3)) 
v2(4)=(v(4)+.25 *abs(v(4)}-v(3))) ; 
v2(3)=v(3), 
else 
v2(3)=w(3); 
v2(4)=v(4); 
end 
end 
aXxis( V2); 
grid on 
hh=slegend((COMPUTED DAMPING ERROR',, TRUE DAMPING ERROR’); 
axes(hh),; 
if pswitch=='y' 
prtfig(fignum) 
delete(h) 
end 
clear DC true_damping true_mass true_stiffness 


YLA%WAVVLVWLYVALVVVVVUVNUVVVYOVAVVVWLVAWVAVWLAWWWWWWVVWVVVWVV%WVV0%Y% 
%%VM% 


load H_EXP 

load CORRHA 

%for mid_index=1:min(cset_size,multplot) 

h=figure(h+1), 

fignum=fignum+] ; 

plot(w/(2*pi),log 1 0(abs(H_EXP(ndx3d([1 aset_size*(aset_size+1)/2 

numpoints], 1 ,skyred(cset_rel(mid_index),cset_rel(mid_index)),"")))),'r--..... 
w/(2* pi), log 10(abs(CORRHA(ndx3d({aset_size aset_size numpoints],cset_rel(mid_index),cset_rel(mid_index),"')))),b-.’) 
vlabel(("H(,int2str( cset(mid_index)),’,,mt2str(cset(mid_index)),') un/Ibf (Log10 of)'}) 
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xlabel((Omega (Hz)'’) 
if titles | 
title(Experimental FRF vs DZ Corrected FRF') 
end 
end on 
V=axIS; 
hh=slegend(EXPERIMENTAL FRF’,'CORRECTED FRF'), 
axes(hh) ; 
if pswitch=="'y' 
prtfig(fignum) 
delete(h) 
end 


clear CORRHA 

load CORRHAD 

h=figure(h+1); 

fignum=fignum+]; 

plot(w/(2*p1), log! 0(abs(H_EXP(ndx3d([1 aset_size*(aset_size+] )/2 

numpoints],1,skyred(cset_rel(mid_index),cset_rel(mid_index)),"')))),'r—",... 
W(1:numpoints-2)/(2*pi),log 10(abs(CORRHAD(ndx3d([aset_size aset_size (numpouits- 

2)],cset_rel(mid_index),cset_rel(mid_index),"’)))),'b-.’) 

ylabel({'H(’,int2str(cset(mid_index)),',",int2str(cset(mid_index)),") in/Ibf (Log10 of)']) 

xlabel(‘Omega (Hz)') 

if titles 
title(Experumental FRF vs DK/DM/DC Corrected FRF') 

end 

erid on 

v=anis; 

hh=slegend(EXPERIMENTAL FRF',,;CORRECTED ANALYTIC FRF’), 

axes(hh); 

if pswitch=='y' 
prtfig(fignum) 
delete(h) 

end 

clear CORRHAD 


save plotnum fignum 
clear 
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MM%MYVSSTCONF FILE FOR A SPATIALLY COMPLETE BEAM 
beanimndl ;%beammdl is a cantilevered 20 element beam model 
%suundl ;%simmdl 1s a masses and springs model 

pswitch=""" —_ ;%%do we print? 

titles=1 ;%%display titles 


meters=1 ;Youse progress meters 
whitebg(‘black’) ;Yoswitch to black figure background 
close 

lowmode=1; 


highmode=1 0; 

static=1;%%%%%N%%%%*H%%V%%reduction method 0O=Guyan, 1=IRS,2=Extraction 

%%define A set & O set%M%%%H%A%%MM%MM%MMMMNMMMMMMMMNMNMMMMMMM%OM% 
aset=1 :numdof: 


%oset=2:2:numdof; %%A SPATIALLY INCOMPLETE BEAM%%%%%%%%%% 
%oset=sort(oset); %%%%%%%%MMUMMUMMMMMMNMUMMMNMMNMUMNMUMNMNMNMNNMUMNMM% 
oset=[]; %%%%A SPATIALLY COMPLETE BEAM%%%%%%%%%%M%%%M%% 

save sstconf 


sst 
plotsst 
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WHY%YSSTCONF FILE FOR A SPATIALLY INCOMPLETE BEAM 
beanundl ‘Ybeanundl is a cantilevered 20 element beam model 
%simundl ;Y%simmdl is a masses and springs model 

pswitch='n' ;%%do we print? 

titles=1 ;%Y%display titles 


meters=1 ;Youse progress meters 
whiutebg(‘black") ;Yeswitch to black figure background 
close 

lowmode=1; 


highmode=10; 


static=1;%%%%%%%M%VM%M%NM%%Y%reduction method 0O=Guyan, l=IRS,2=Extraction 

%%define A set & O set%ww%%MM%%M%%%%%HM%%MMMMMM%%%%MNMNM%MNM%%M%%M% 
aset=] :numdof; 

oset=2:2:nuindof, %%A SPATIALLY INCOMPLETE BEAM%%%%%%%%%% 
oset=sort(oset); %%%%%%%%%MH%%%%MMMMMMMMM%MNMMVMU%%MNM%%V%ON% 


save sstconf 
sst 
plotsst 


beammad] -Ybeanimndl is a cantilevered 20 element beam model 
%simmd] ;Y%sunmdl is a masses and springs model 

pswitch='n' ;%%do we print? 

titles=1 ;%Y%display titles 


meters=] ;Youse progress meters 
whitebg(‘black’) ;Yoswitch to black figure background 
close 

lowmode=1; 


highmode=10; 


%Ydefine A set & O set%>%%%%%%%M%%%MNMM%M%MMH%M%MUMNM%MUHNMUMM%NNMMN%N% 
aset= 1 :numdof; 
oset=2:2:nuindof, %%A SPATIALLY INCOMPLETE BEAM%%%%%%%%%% 
oset=sort(oset); H%%%%%N%%%M%%%M%MMMMMMM%MNMMMNMMMNMMMMMMNMNMNM%NN% 


static=2; 


save sstconf 
dofig4 3 
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WWMWWW“MVVHVM“VUVUVMUWUVUVUMVWVUVWYVVWMUVM%M%MMVWVNMUMAMNMVVMWMYUVVVUNNMWVWNWWPVYUWWYWYVWWVY% 
WUWUMVWM%wWVM“MVUN% 

%%% COMPARE SINGLE MODE MATRIX AND INTEGRAL %*%%%%KUUM“VKHYUYwWNYYWYMVVYVUYWWAHVVY 
%%% SOLUTIONS UNDER VARIOUS CONDITIONS  %%%%W%%%%MUWMKUYUwHUVUMNMVVYnHYVNWVUVWVYV% 
UWUVW%V%%VVAV%NM%MUVVVNMNVNMUMMNMUVNMUWNMMMNMMMMNMMMVMNMMMNMUMNMUWMNMNVUVV“VUVYUVUYUWVUYUVVNVVV% 
BIL%%NM%MNMMM%M%O% 

clear 

closeall 

load INT 

whitebg(‘white’) 

clear L Z_ ANAL c_anal c_exp conn h_exp h_anal_red k_anal k_exp 

clear kexta kexto m_anal m_exp mexta mextO temp_l_diags 

clear true_damping true_mass true_stif{ness z_anal z_anal_red z_exp 

clear L Z_ANAL c_anal c_exp conn h_exp h_anal_red k_anal 

clear kexta kexto m_anal mexta inext0 temp_l_diags 


load H_EXP 

load exp 

cset=[5 79 11 13 15] 

cset_rel=[3 45 67 3] 

cset_size=length(cset) 

muid_index=round(length(cset)/2); 

fineness=200 

ODZ1=[]: 

temp1=[]; 

odtreq=omegax(1); 

center_freq=odfreq/(2*p1) 

w1=odfreq; 

odlength=2* pt; 

odivisions=2; 

lowerfreq=odfreq-15*odlength 

upperfreq=odfreq+25*odlength 

ow]1=odfreq-.5* odlength:odlength/odivisions:odfreqt+.5* odlength; 

partition=ow1/(2* pi) 

df=(odlength/odivisions)/(2*p1) 

for count=1:length(ow1) 
Z_anal_red=kstat-owl(count)*2*mstat; 
h_anal_red=inv(z_anal_red), 
Z_exp=k_exptj*owl(count)*c_exp-owl (count)*2*m_exp, 
h_exp=inv(Z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcec=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hace-hxcc)* inv(hacc))-hacc); 
ODZ1=[ODZ1, dz); 
temp1=[temp1; eye(cset_size) -owl(count)*2*eye(cset_size) +j*ow](count)*eye(cset_size)}; 


end 

“%theleftdz=ODZ1 

%thebigone=temp 1 

DKDMDC1=temp1\ODZ1; 
ODK1=DKDMDC1(1:cset_size,:) 
ODM1=DKDMDC 1 (cset_size+1:2*cset_size,:) 
ODC1=DKDMDC1(2*cset_size+1:3*cset_size,:) 
save ODM1 ODM! DKDMDC1 ODZ1 owl wl 
save ODK1 ODK1 

save ODC1 ODC1 


Yo HH 
load stat; 
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odcorrha={]; 
w=lowerfreq:(upperfreq-lowerfreq)/funeness:uppertfreq; 
numpoints=length(w),; 
kcorrected=kstat; 
65 imcorrected=mstat; 
ccorrected=cstat; 
kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel }ODK 1; 
meorrected(cset_rel,cset_rel)=incorrected(cset_rel,cset_rel}+ODM1; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel}+ODC1; 
70 for 1=1:numnpoits 
tempza=kcorrected+j*w(i)* ccorrected-w(1)*2 *mcorrected; 
tempha=inv(tempza); 
odcorrha=[odcorrha tempha(:)]; 
end 
75 
uncorrha=[]; 
for i=1:numpoints 
tempza=kstat+j*w(1)*cstat-w(1)*2 *imstat; 
tempha=inv(tempza), 
80 uncormha=[uncorrha tempha(:)]; 
end 


if meters 
waitbar_handle=waitbar(0,'Computing Experimental FRF"); 
85 else 
disp(Getting experimental FRF") 
end 
H_EXP=[]; 
for count=1:length(w) 
90 if meters 
waitbar(count/length(w)); 
end 
z_exp=k_expt+j*w(count)*c_exp-w(count)*2*m_exp; 
h_exp=inv(z_exp); 
95 h_exp=h_exp(aset,aset); 
skyindex=1; 
for index1=] :aset_size 
for index2=index]:aset_size 
ted_holder(skyindex)=h_exp(index 1 ,index2); 
100 skyindex=skyindex+ 1; 
end 
end 
H_EXP=[H_EXP red_holder’ 
end 
105 if meters 
close{waitbar_handle) 
end 
clear tempha tempza kcorrected mcorrected ccorrected 


110 figure( 1); 
plot(w/(2*pi),log10(abs(H_EXP(ndx3d([1 aset_size*(aset_size+1)/2 numpounts].... 


1,skyred(cset_rel(mid_index),cset_rel(imid_index)),"")))),'-’.... 
w/(2*pi),log 10(abs(odcorrha(ndx3d([aset_size aset_size nuinpoints],cset_rel(mid_imdex),cset_rel(mid_index),'')))),'’g—'.... 


w/(2*p1),log 10(abs(uncorrha(ndx3d((aset_size aset_size numpoints]},cset_rel(mid_index),cset_rel(mid_index),’:')))),'b:’) 
ylabel({'H(,mt2str(cset(mid_index)),',’,int2str(cset(mid_index)),") m/lbf (log10 of)']) 


120 xlabel('Omega (Hz)’) 
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“%title(‘Single Mode Corrected FRF (Mode 1) ') 
grid 
hh=slegend(‘Experimental FRF', ‘Mode 1 Corrected FRF', ‘Uncorrected Analytic FRE"); 


axes(hh) 
print -dinfile fig5_1 


AWRwWVVVVUwVMVVMUVUWVMUVMUUVUWMVUUUwUMUUVUVUVUVUVUVUVUUUVHVUVUVUVUVUVUVAV*VM“VUVUVUwVUMUAVUUAUVUVUUHUV% 
YHYHDEOOMPOBRVDOEL into DM, DK, & DC RRwR%RwW%*X*%VZ%**R%MU“MZwVUVUV%UUVUUwMUMUVMLUVMVwV“VUVAMVMUUMUMUVHUV%VN% 
%*wY% USING MATRIX TECHNIQUE %HW%UW%%%%M%%M%U%%UUVVUVWUUMVUNUVVU%VMVUUAUM%VUV%% 
RWVAV/HV“LVMVHAHVVUwNY*WMVUVV“HVVVMVHVVWwVU%wWVUUWVYVUVVVYUVVU“VVUVVVUVU*VUVUHUYUHUUVUUMUVUWN% 
QDPAHN%%%NMNMUN 

temp 1={]; 

odlength=2 *pi; 

owl=oniegax(1 )-5*odlength:(omegax(2)-omegax(1)+10* odlength)/14:omegax(2)+5*odlength; 


for count=1:length(owl ) 
Z_anal_red=kstat-owl(count)*2*mstat; 
h_anal_red=inv(z_anal_red), 
Z_exp=k_exptj*owl(count)*c_exp-owl (count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc )*(hacc-hxec)*inv(hacc))-hacc); 
ODZ1=[ODZ1; dz}; 
temp1=[temp1; eye(cset_size) -owl(count)*2*eye(cset_size) +j*ow](count)*eye(cset_size)}; 


end 


DKDMDC1=temp1\ODZ1, 

ODK1=DKDMDC 1(1:cset_size,:) ; 
ODM1=DKDMDC 1 (cset_size+1:2*cset_size,:) ; 
ODC1=DKDMDC1(2*cset_sizet1:3*cset_size,:); 
save ODM1 ODM1 DKDMDC1 ODZ1 owl wl 
save ODK1 ODK1 

save ODC1 ODC1 


UWNWWVAWHVwWV%MNMM%UWMN%MMWMHWNMNMUWMMNNMVNHNVWAVMVW%MUWM%%%MMMAV%MVAVWNMMV%N% 
PORIKMIY% MMMM 
odcorrha=[]; 
w=omegax(1)-15*pi:(omegax(2 }-omegax(1 +35 *pi)/fineness:omegax(2 }+20*p1; 
numpoints=length(w); 
kcorrected=kstat; 
incorrected=mistat; 
ccorrected=cstat; 
kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel}}+}ODK 1; 
mcorrected(cset_rel,cset_re] )}=mcorrected(cset_rel,cset_rel}}ODM1,; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel}}ODC1; 
for i=1:numpoints 
tempza=kcorrected+j*w(i)*ccorrected-w(i)*2*mcorrected; 
tempha=inv(tempza); 
odcorrha=[odcorrha tempha(:)]; 
end 


uncormha={]; 

for i=] :numpoints 
tempza=kstat+j * w(i)*cstat-w(i)*2*mstat; 
tempha=inv(tempza): 
uncorrha=[uncorrha tempha(:)}, 


108 


185 


190 


195 


200 


205 


210 


220 


geo 


230 


235 


240 


CHAPS.M 


end 


if meters 
waitbar_handle=waitbar(0,'Computing Experimental FRI"); 
else 
disp(‘Getting experimental FRI") 
end 
H_EXP={]; 
for count=1:length(w) 
if meters 
waitbar(count/length(w)); 
end 
Z_exp=k_exptj*w(count)*c_exp-w(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex=1; 
for index1=1:aset_size 
for index2=index | :aset_size 
red_holder(skyindex)=h_exp(index 1 ,ndex2); 
skyindex=skyindex+1; 
end 
end 
H_EXP=[H_EXP red_holder’}; 
end 
if meters 
close(waitbar_handle) 
end 
clear tempha tempza kKcorrected mcorrected ccorrected 


figure(2); 
plot(w/(2*pi),log10(abs(H_EXP(ndx3d([1 aset_size*(aset_size+])/2 numpoints].... 


1,skyred(cset_rel(mid_index),cset_rel(imid_index)),"’)))),'r~",... 
w/(2*pi),log10(abs(odcorrha(ndx3d([aset_size aset_size numpoints],cset_rel(mid_index),cset_rel(mid_index),"')))),'g—.... 


wh(2*pi),log10(abs(uncorrha(ndx3d({aset_size aset_size numpoints],cset_rel(mid_index),cset_rel(mid_index),"’)))),b") 
ylabel(("H(,int2str(cset(mid_index)),",",int2str(cset(mid_index)),') in/Ibf (log10 of)']) 


xlabel('Omega (Hz)’) 

Ytitle(‘Single Mode Corrected FRF (Mode 1) ') 

grid 

hh=slegend(Experimental FRF’, ‘Single Mode Corrected FRE’, ‘Uncorrected Analytic FRF’); 


axes(hh) 
print -dinfile fig5 2 


WHVAWV/HVVHAVVVVAVUVAVLVAVMHVUWHVUUVAHUVUVVVVUVVAVUUWVUUUAHUVVVUVUVVVVYHWYVY 
YWVV OPasenMmreDsA nto ODM, ODK, & ODC %%*%%%w%M%MNNMMUMUMUMUUUUUUMUVUUNMUVVVVO% 
R~WR~VMNVUWUVVLLU“MMUUwUUMUANMUMUVUMVAUVAHUMUNMUWUVUWVAVUHUVVU“VWUUUVUVVAVUUVUUUVVUNNWVYWUVUWYV 
FABER Sas So iG Loisd Pi eq-lowerfreq)/fineness:upper freq; 
lenctrl=[(25 10 1); 
odivisions=3; 
for ind=1:length(lenctrl) 

odlength=lenctrl(ind)*2*p1; 

ODZ=(); 

temp=(]; 

ow=odfreq-.5*odlength:odlength/odivisions:odfreq+.5*odlength; 

for count=] -length(ow) 

z_anal_red=kstat-ow(count)*2*mstat; 
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h_anal_red=inv(z_anal_red), 
Z_exp=k_exptj *ow(count)*c_exp-ow(count)*2*m_ exp; 
h_exp=inv(z_ exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel), 
hxcc=h_exp(cset_rel,cset_rel); 
=inv(inv(inv(hacc)*(hace-hxcc)* inv(hacc))-hacc); 
ODZ=[ODZ; dz]; 
temp=[temp; eye(cset_size) -ow(count)*2*eye(cset_size) +)*ow(count)*eye(cset_size)]; 


end 


DKDMDC=temp\ODZ; 
ODK=DKDMDC(1:cset_size,:) ; 
ODM=DKDMDC(cset_sizet1:2*cset_size,:) ; 
ODC=DKDMDC(2* cset_size+1:3*cset_size,:); 
clear hace hxcc h_exp z_exp z_anal_red h_anal_red dz 
save ODM ODM DKDMDC ODZ ow 
save ODK ODK 
save ODC ODC 
WUWwWNYwAVwVMAwHUVM’VNMwM*MwMM*wMMnHUwMUVUMNMMMMUMMNNMAMnMMMNMMNMUMMUNM*MAUMUANMU%NMUM%M%UMNNMNUN%N% 
PAV AEP %%%% 
waitbar_handle=waitbar(0,'Computing Experimental FRF’), 
else 
disp('Getting experimental FRF") 
end 
H_EXP={]; 
for count=1 :length(plotw) 
if meters 
waitbar(count/length(plotw)); 
end 
Z_exp=k_expty*plotw(count)*c_exp-plotw(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex=1; 
for index1=1:aset_size 
for index2=index] :aset_size 
red_holder(skyindex)=h_exp(index1 ,index2); 
skyindex=skyindex+ 1; 
end 
end 
H_EXP=(H_EXP red_holder'}; 
end 
if meters 
close(waitbar_handle) 
end 


odcorrha=[]; 
kcorrected=kstat; 
mcorrected=mstat; 
ccorrected=cstat; 
kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel}tODK; 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel }}ODM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel HODC; 
for 1=1:length(plotw) 
tempza=kcorrected+ *plotw(1)*ccorrected-plotw(1)*2*mcorrected; 
tempha=inv(tempza); 
odcorrha=[odcorrha tempha(:)]; 
end 
if ind==1 
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plotl=odcorrha; 
elseif ind== 
plot2=odcorrha; 
elseif ind==3 
plot3=odcorrha; 
else 
plot4=odcorrha; 
end 
end 


figure(3), 
plot(plotw/(2 *pi),log10(abs(H_EXP (ndx3d({1 aset_size*(aset_size+1)/2 length(plotw)].... 


1 skyred(cset_rel(mid_index),cset_rel(muid_index)),":’)))),'r-\... 
plotw/(2* pi), log10(abs(plot1(ndx3d([aset_size aset_size length(plotw)],cset_rel(imid_index),cset_rel(imid_index),'')))),'¢~ 


plotw/(2*pi), log10(abs(plot2(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),":’)))),'b’... 
plotw/(2*pi), log10(abs(plot3(ndx3d([aset_size aset_size length(plotw)],cset_rel(mid_index),cset_rel(mid_index),'')))),‘k-.’) 


ylabel({'H(‘,int2str(cset(mid_index)),',’,int2str(cset(mid_index)),') in/lbf (log10 of)']) 


xlabel(‘Omega (Hz)’) 
gnd on 
%title"EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 5 SOLUTIONS’) 
hh=legend(‘Exp FRF*, ... 
{num2str(lenctrl(1)),' Hz Bandwidth Corrected FRF', ... 
[num2str(lenctrl(2)),' Hz Bandwidth Corrected FRF'J, ... 
{num2str(lenctrl(3)),' Hz Bandwidth Corrected FRF']) ; 
axes(hh) 
print -dinfile fig5_3 


odlength=1*2*p1; 
countctrl=[3 10 50 200] ; 
for ind=1:length(countctrl) 
odivisions=countctrl(ind)-1; 
ODZ=(]; 
temp=[]; 
ow=odfreq-.5* odlength:odlength/odivisions:odfreqt+.5*odlength; 
for count=1:length(ow) 
z_anal_red=kstat-ow(count)*2 *mstat; 
h_anal_red=inv(z_anal_red), 
Z_exp=k_exptj*ow(count)*c_exp-ow(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hace-hxcec)*inv(hacc) -hacc); 
ODZ=[ODZ,; dz]; 
temp=[temp; eye(cset_size) -ow(count)*2*eye(cset_size) +j}*ow(count)*eye(cset_size)]; 


end 


DKDMDC=temp\ODZ; 
ODK=DKDMDC(l:cset_size,:) ; 
ODM=DKDMDCcset_size+1:2*cset_size,:) ; 
ODC=DKDMDC(2*cset_size+1:3*cset_size,:); 

clear hace hxcc h_exp z_exp z_anal_red h_anal_red dz 
save ODM ODM DKDMDC ODZ ow 

save ODK ODK 


111 


365 


370 


375 


380 


385 


390 


395 


405 


410 


415 


420 


CHAPS.M 


save ODC ODC 
~P~YYRL/HVAVVVHL“HVAUVANMM%W%%M%MNMMMMMMNUMUMMWMNMMMMVWWVWVNWMVNMU%VMMMNMANMUMUVUMMUVU%UN% 
wV% 
odcorrha=[]; 
kcorrected=kstat; 
incorrected=mstat; 
ccorrected=cstat; 
Kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel HODK; 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_re] +} ODM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel #ODC; 
for i=1:length(plotw) 
tempza=kcorrected+)* plotw(1)*ccorrected-plotw(i)*2 *mcorrected; 
tempha=inw(tempza); 
odcorrha=[odcorrha tempha(:)]; 
end 
if ind==1] 
plot 1=odcorrha; 
elseif ind== 
plot2=odcorrha; 
elseif ind==3 
plot3=odcorrha; 
else 
plot4=odcorrha; 
end 
end 


figure(4); 
plot(plotw/(2*p1),log10(abs(H_EXP(ndx3d({1 aset_size*(aset_size+1)/2 length(plotw)].... 


1 ,skyred(cset_rel(mid_index),cset_rel(mid_index)),"’)))),'r-... 
plotw/(2*p1), log] 0(abs(plot 1(ndx3d({aset_size aset_size length(plotw)],cset_iel(mid_index),cset_rel(mid_index),"')))),’g— 


plotw/(2*pi), log10(abs(plot2(ndx3d((aset_size aset_size 
length (plotw)],cset_rel(imid_index),cset_rel(mid_mdex),"’)))),’b’’... 
plotw/(2*pi), log 10(abs(plot3(ndx3d({aset_size aset_size length(plotw)],cset_rel(mid_index),cset_rel(mid_index),"’)))),‘k- 


| plotw/(2*pi), log 10(abs(plot4(ndx3d({aset_size aset_size length(plotw)],cset_rel(mid_index),cset_rel(mid_index),'"’)))),m°) 
ylabel({'H(,int2str(cset(nud_index)),',\int2str(cset(nud_index)),") in/lbf(10og10 of)']) 


xlabel(‘Omega (Hz)') 
grid on 
%Mtitle(EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 1 SOLUTIONS’) 
hh=legend(‘Exp FRF*, ... 
{num2str(countctrl(1)),' Pts Corrected FRF’], ... 
{num2str(countctrl(2)),' Pts Corrected FRI", ... 
{nun2str(countctrl(3)),’ Pts Corrected FRF’'J, ... 
{num2str(countctrl(4)),’ Pts Corrected FRF']); 
axes(hh) 
print -dinfile figS 4 


MVWVVVVVHVMHKHM%M%%M%MMMMMMMMMMM%MW%%M%MM%%MMMMMMMMMMMMMMMMMMMNMMMM%MUMM% 
WV/HDELOMPORMYDZ into DM, DK, & DC %%%%%%%%%NMMNMMMMMMMWNMMMUMNM’MNMNMMUMUMUNMNHUV%N% 
%%Y% USING INTEGRAL TECHNIQUE %%%%%%%%%%%M%%%MMM%M%MMMMMMMMMUMMUMMNMMMNM%N% 
~AVALLHhHHAMANVVKA“wVNWM’MKVUUMMUAMMMMNMM%MMMMMMMMMMMNMMMMMMMMMUMMNUMNMUMUMM% 
Ye PoDplo%%%%% 

Y%plotw=lowerfreq:dw-upperfreq; 

lenctrl=[25 10 1]; 

odivisions=3; 

for ind=1:length(lenctr]) 
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odlength=lenctrl(ind)*2 *p1; 

owl=odfreq-.5 *odlength:odlength/(odivisions):odfreqt. 5 *odlength; 

Wl=owl; 

%YYocompute integrals YWUWVUVMVUKH“VUMMUMUVHUVAVHMUVVWVVVUWVVWUVVWVUWVVVVVV%% 

Yweight=ones(size(ow] ))*(omegax(1 )); 

weight=ow]; 

DZ_R=[]; 

DZ_I=(]: 

for count=1:length(ow1) 
z_anal_red=kstat/(j*ow1(count)}+j*ow1(count)*mstat; 
h_anal_red=inv(z_anal_red), 
z_exp=k_exp/(j*ow](count))}+c_exp/(j*ow1(count))+j*owl(count)*m_exp; 
h_exp=inv(z_exp), 
h_exp=h_exp(aset,aset); 
hace=h_anal_red(cset_rel,cset_rel); 
hxec=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hacc-hxcc)*inv(hacc))}-hacc); 
dz_r=real(dz); 
dz_i=1mag(dz); 
DZ_R=(DZ_R dz_r(:)}; 
DZ_I=[DZ_I dz_i(:)]; 

end 


(INTK INTM INTC}=intsub(DZ_LDZ_R,ow1,W1,1./4 *ow] ),cset_size); 
odcorrha={]J; 
load stat 
kcorrected=kstat; 
mcorrected=mstat; 
ccorrected=cstat; 
kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel +INTK; 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel HINTM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel H INTC; 
for 1=1:length(plotw) 
tempza=kcorrected+j*plotw(1)*ccorrected-plot w(1)*2*mcorrected; 
tempha=inv(tempza), 
odcorrha=[odcorrha tempha(:)}; 
end | 
if ind== 
intplot1=odcorrha; 
save INTM] INTM 
save INTK] INTK 
save INTC] INTC 
elseif ind== 
intplot2=odcorrha;, 
save INTM2 INTM 
save INTK2 INTK 
save INTC2 INTC 
elseif ind== 
intplot3=odcorrha; 
save INTM3 INTM 
save INTK3 INTK 
save INTC3 INTC 
else 
intplot4=odcorrha; 
save INTM4 INTM 
save INTK4 INTK 
save INTC4 INTC 
end 
end 
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if meters ' 
waitbar_handle=waitbar(0,'Computing Experimental FRF’); 
else 
disp(‘Getting experimental FRF') 
end 
H_EXP={]; 
for count=1:leugth(plotw) 
if meters 
waitbar(count/length(plotw)); 
end 
Z_exp=k_exptj*plotw(count)*c_exp-plotw(count)*2*m_exp; 
h_exp=mv(z_ exp), 
h_exp=h_exp(aset,aset), 
skyindex=1, 
for index1=1:aset_size 
for index2=index] :aset_size 
red_holder(skyindex)=h_exp(index 1 ,index2),; 
skyindex=skyindex+1; 
end 
end 
H_EXP=[H_EXP red_holder’}; 
end 
if meters 
close(waitbar_handle) 
end 


figure(5); 
plot(plotw/(2* pi),log 1 0(abs(H_EXP(ndx3d([1 aset_size*(aset_sizet+1)/2 length(plotw)].... 


1 ,skyred(cset_rel(muid_index),cset_rel(mid_index)),"')))),'r-.... 

plotw/(2*pi), log10(abs(intplot 1(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),":')))),'g—'.... 

plotw/(2*pi), log10(abs(intplot2(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),"")))),'b".... 

plotw/(2*pi), log10(abs(intplot3(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),"')))),k-.’) 
ylabel((H(‘,int2str(cset(imid_index)),',",int2str(cset(mid_uidex)),") in/lbf (1og10 of)']) 


xlabel(‘(Omega (Hz)') 
grid on 
%title"EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 1 INT SOLUTIONS’) 
hh=legend(‘Exp FRF', ... 
[num2str(lenctrl(1)),' Hz Bandwidth Corrected FRF').... 
[num2str(lenctrl(2)),' Hz Bandwidth Corrected FRF'].... 
[{num2str(lenctrl(3)),' Hz Bandwidth Corrected FRF'}); 
axes(hh) 
print -dinfile fig5 5 


odlength=2 *pi; 
countctrl=[3 10 50 200}; 
for ind=1 :length(countctrl) 
odivisions=countctrl(ind); 
ow]=odfreq-.5* odlength:odlength/(odivisions-1):odfreq+.5*odlength; 
W1=owl; 
%Y%compute integrals %%%%%%%%H%%%%%M%%%%WHM%MMM%%%%%M%%%%%MM%%%%%%M%M%%V0% 
“%weight=ones(size(ow 1 ))*(omegax(1)); 
weight=ow1; 
if ind == 
owl 
df=(odlength/(odivisions-1))/(2* pi) 
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end 

%% compute integrals M%wW%wR%wR%W%wYW%°%V%V%%UUUUUURUU*VVwV*Vw“wVwVUVUHUwWwMUwUV%UUwMUVUVUVUYUMVMUWM% 
DZ_R={]; 

DZ_I={]; 


for count=1-:length(ow1) 
z_anal_red=kstat/(j **ow1(count)}+j*ow1(count)*mstat; 
h_anal_red=inv(z_anal_red); 
z_exp=k_exp/(j*owl(count)}+c_exp/(j*owl (count) )}+j *owl(count)*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc )*(hacc-hxcc)*inw(hacc ))-hacc); 
dz_r=real(dz), 
dz_i=1mag(dz); 
DZ_R=[DZ_R dz_x(:)]; 
DZ_I=[DZ_I dz_i(:)]; 
end 
if ind == 
DZ_R; 
DZ_I; 
for 1=1:length(owl ) 
dztemp1=DZ_R(ndx3d({1 cset_size*(cset_size+1)/2 length(ow1)],1,1:cset_size*(cset_sizet+1 )/2,1)); 


dztemp2=DZ_I(ndx3d([1 cset_size*(cset_size+1)/2 length(ow1)],1,]:cset_size*(cset_size+])/2,i)); 


skyindex=1; 
for index1=1:cset_size 
for index2=index 1 :cset_size 
dz1(index 1 ,index2)=dztemp1(skyindex), 
dz1(index2,index1)=dztemp1(skyindex), 
dz2(index 1 ,index2)=dztemp2(skyindex), 
dz2(index2,index1)=dztemp2(skyindex); 
skyindex=skyindex+ ; 
end 
end 
dzl 
dz2 
end 


end 
[INTK INTM INTC]}=intsub(DZ_LDZ_R,ow1,W1,1./G *ow1),cset_size); 
if ind == 

(INTK INTM INTC]=indsub(DZ_LDZ_R,ow1,W1,1./G*ow1),cset_size); 
INT™M 


INTK 
INTC 
end 


YRS lt pP/0% 7 %% UMM %%%%VVWUW%%UMMVMUVVAVWMUAVUUMUMVMMUVMVUVMAMUVMMUMUVNWMVUMVWUVN% 
YALarhstat 
kcorrected=kstat; 
mcorrected=mstat; 
ccorrected=cstat; 
kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel +INTK, 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel + INTM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel +INTC; 
for i=1:length(plotw) 
tempza=kcorrected+j*plotw(1)*ccorrected-plotw(1)*2 *mcorrected; 
tempha=inv(temipZa),; 
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odcorrha={odcorrha tempha(:)]; 

end 

if ind== 
mtplot1=odcorrha; 
save INTM1 INTM 
save INTK1 INTK 
save INTC! INTC 

elseif ind==2 
intplot2=odcorrha; 
save INTM2 INTM 
save INTK2 INTK 
save INTC2 INTC 

elseif ind== 
intplot3=odcorrha; 
save INTM3 INTM 
save INTK3 INTK 
save INTC3 INTC 

else 
intplot4=odcorrha; 
save INTM4 INTM 
save INTK4 INTK 
save INTC4 INTC 

end 

end - 


if meters 
waitbar_handle=waitbar(0,'Computing Experimental FRE’), 
else 
disp(‘Getting experimental FRF’) 
end 
H_EXP=[]; 
for count=1:length(plotw) 
if meters 
waitbar(count/length(plotw)); 
end 
z_exp=k_exptj*plotw(count)*c_exp-plotw(count)*2*m_exp; 
h_exp=imv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex=1, 
for index1=1:aset_size 
for index2=index]:aset_size 
ted_holder(skyindex)=h_exp(index] ,index2); 
skyindex=skyindex+1, 
end 
end 
H_EXP=[{H_EXP red_holder’}; 
end 
if meters 
close(waitbar_handle) 
end 


figure(6); 


plot(plotw/(2*pi),Jog10(abs(H_EXP(ndx3d({1 aset_size*(aset_size+1)/2 length(plotw)].... 


1,skyred(cset_rel(mid_index),cset_rel(mid_index)),"')))),'r-".... 

plotw/(2*pi), log10(abs(intplotl(ndx3d({aset_size aset_size 
length(plotw)]},cset_rel(mid_index),cset_rel(mid_index),""')))),'g—".... 

ploitw/(2*pi), log10(abs(intplot2(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),'.’)))),"b-.... 
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plotw/(2*pi1), log10(abs(intplot3(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),':)))),'k-.‘... 

plotw/(2*p1), log10(abs(intplot4(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),':')))),'m:’) 
ylabel(["H(‘,int2str(cset(mid_index)),',int2str(cset(inid_index)),') in/Ibf (log10 of)']) 


xlabel(‘Omega (Hz)’) 
erid on 
Ytitle"EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 1 INT SOLUTIONS’) 
hh=slegend(‘Exp FRF’, ... 
[num2str(countctrl(1)),' Pts Corrected FRF', ... 
[nun2str(countctrl(2)),' Pts Corrected FRF', ... 
{num2str(countctrl(3)),' Pts Corrected FRF’], ... 
{numn2str(countctr](4)),' Pts Corrected FRF’]); 
axes(hh); 
print -dinfile fig5 6 
figure(7); 
plot(plotw/(2*pi),log10(abs(H_EXP(ndx3d([1 aset_size*(aset_size+1)/2 length(plotw)].... 


1 ,skyred(cset_rel(mid_index),cset_rel(mid_index)),"')))),'r-'.... 
plotw/(2* pi), 1og10(abs(plotl(ndx3d([aset_size aset_size length(plotw)],cset_rel(mid_index),cset_rel(imid_index),':')))),'g— 


1 
>° 


plotw/(2*p1), log10(abs(intplot 1(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(imid_index),‘')))),'b:') 
ylabel((HC,int2str(cset(mid_index)),",",int2str(cset(mid_index)),") uV/lbf (log10 of)']) 


xlabel(Omega (Hz)') 

if titles 

% title(EXPERMENTAL FRF AND CORRECTED FRFS USING MODE | OD&INTGRL SOLTNS’') 

end 

erid on 

hh=legend(‘Exp FRF’, ... 
[num2str(lenctri(1))," Hz Bandwidth MAT Corrected FRF'J, ... 
[num2str(lenctrl(1)),' Hz Bandwidth INT Corrected FRF); 

axes(hh) 

if pswitch=='y' 

print -dedjcolor 
end 


print -dinfile figS_7 


figure(8); 
plot(plotw/(2* pi),log10(abs(H_EXP(ndx3d([1 aset_size*(aset_size+1)/2 length(plotw)].... 


1 ,skyred(cset_rel(mid_index),cset_rel(mud_index)),"")))),'r-"’,... 
plotw/(2*pi), log10(abs(plot2(ndx3d([aset_size aset_size length(plotw)],cset_rel(mid_index),cset_rel(mid_index),"’)))),'g— 


plotw/(2*pi), log10(abs(intplot2(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mid_index),"“")))),’b:') 
ylabel(['H(‘,int2str(cset(mid_index)),',,int2str(cset(mid_index)),") i/Ibf (log10 of)']) 


xlabel((Omega (Hz)') 
if titles 
%title"EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 1 OD&INTGRL SOLTNS') 
end 
gnd on 
hh=legend(Exp FRF', ... 
{num2str(lenctrl(2))," Hz Bandwidth MAT Corrected FRF'J, ... 
{num2str(lenctrl(2)),' Hz Bandwidth INT Corrected FRF']); 
axes(hh) 
if pswitch=='y' 
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print ~dcedjcolor 
end 
print -dinfile fig5 8 


figure(9); 
plot(plotw/(2* p1),log 1 0(abs(H_EXP(ndx3d([1 aset_size*(aset_size+])/2 length(plotw)}.... 


1 ,skyred(cset_rel(mid_index),cset_rel(mid_index)),""')))),'r-’,... 
plotw/(2*pi), log 10(abs(plot3(ndx3d({aset_size aset_size length(plotw)},cset_rel(imid_index),cset_rel(mid_index),‘:')))),'g— 


plotw/(2*pi), log10(abs(intplot3(ndx3d([aset_size aset_size 
length(plotw)],cset_rel(mid_index),cset_rel(mud_index),'')))),"b:') 
ylabel(('H(,int2str(cset(mid_index)),',"int2str(cset(mid_index)),") in/Ibf (log10 of)']) 


xlabel’(Omega (Hz)') 
if titles 
%title( EXPERMENTAL FRF AND CORRECTED FRFS USING MODE 1 OD&INTGRL SOLTNS') 
end 
grid on 
hh=legend(‘Exp FRF’, ... 
{num2str(lenctrl(3)),' Hz Bandwidth MAT Corrected FRF'}, ... 
(num2str(lenctrl(3)),' Hz Bandwidth INT Corrected FRF')); 
axes(hh) 
if pswitch=='y’ 
print -dedjcolor 
end 
print -dinfile fig5_9 
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clear 

odivisions=3 

startloop=1 

skiploop=1 

endloop=4 

cd incomp1 

save odconf odivisions startloop skiploop endloop 
chap6_1 

chap6 2 

cd.. 

clear 

odivisions=| 

startloop=3 

skiploop=1 

endloop=4 

ed comp! 

save odconf odivisions startloop skiploop endloop 


chap6_] 
chap6_3 
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WVWYCHAPH 1 MWMMRwVVMUVVMUMUVMMUVMUVUMUVVUVVHUVVUwY“VUVUVMMMUHM“wVUUVMUMUVAVHVMVHANM“AVMVHAVMUUHNVMUM%MUV%N% 
AYMHAV“V“HVN% 

%w%%Decompose DZ into DM, DK, & DCRV*%*%%%%NM%wMRwR“NMWwMUwMUAVMAUMUMUMAUMNMNMMMUMMNMNMMNMUYRHUVHVMN 
%%wY%using matrix formulationw%w%w%*%%V%~%VV%VnAVHMUwMUUV“VUwMVVMUwVMUVwYUMANVMUVUAAM“MAVVV“LAVUVUAVVHNN% 
WRVVMUM“M’VVWANVVV-“VUUMVMALVHKLANH“VMNMUVUVwMMUMNMUAUVMVHVMVUM“LUVVHVU>MUHUAVLAV“M“VUYANUVVVN% 
YRVW“KAAvHAwMVM%wM%% 

cle 

clear 

closeaill 

load INT 

%if pswitch==' 

whitebg(‘white') 

close 

Y%end 

h=0; 

fignum=1; 

use_antires=1; 

titles=0, 

mid_index=round(length(cset)/2); 


clear L Z_ ANAL c_exp conn h_anal_red 

clear kexto m_exp mexto temp_|_diags 

clear true_damping true_mass true_stiffness z_anal z_anal_red z_exp : 
clear LZ ANAL c_exp conn h_anal_red 

clear kexta kexto mext0 temp_1_diags 


ARWYRWH“HV“K“HVAwMNM*wM“HV“V“U“V“V“%VFORCE THE CSETHR*X*MM%VM*w“VUHU%VWwUMN%NMVO% 
AHRWVHV“WV“VK“VwV*“WYVSPATIALLY INCOMPLETE%%%%UN*VWMVUUwV“MMUHUVV% 
if complete 
WYARVWWMwMKHVVWwUMUAVAKHUV“LVVUUMMMUMUwUMUVUWMVMVMUVMUMMUUUMMUMVMUAVUVUVULUMUMAHVUVVN% 
AYMRwHUV“wMVM%VMVMUwUMUVUH“VVSPATIALLY COMPLETE%%%%%%M%%M%M%MMMMMNN% 
cset=[4 56789 10 1LIJMMMMM>nM*HMNMNNMUMRHU“HUYwVVHUWNY%M% 
cset_rel=csetrMnw%r%UVV%VUHUVMUVUMUV“VUVAHVHVLLUVWUWUUWUVUUMVWUVKVKVWWYY 
cset_size=length(cset )%*n%Z“ZA%KUUWMUMUVUVUMvHUV“LVU“YUVYUMAHVUNNYM 
ARwYwVKLVVHVAwVVM“H“VV“VKw“H“VUVKVHVL“MAMVVUwUUMUV“V“VUVM*wV“VKVMVUVVMMUMUVMVKwHVNVNM% 
else 
RWAwVK“K“K“V*wM*V%“SPATIALLY INCOMPLETE%%%%%%%%ULUMM%UMUHUUM%N% 
cset=[1 3 5 7911 13 ISJ%URK%M%V%W%UWHUwVMVMUUMUVUWUVUVUVU%UUVUVN% 
cset_rel=[1 2 3 45 6 7 8]INRKRwKKKKKKK“KUVUWU~%UNV%V%% 
cset_size=length(cset )%%%HAVKKHKHUVVUwAw“AVVHVWVYHVVYON 
end 
AWRHYVVVMAHKAHK“AHUwHUwMVVVVHVU>wUVwUUVUVVVVwMUVL“VAVL“YVUMAUAMMUHVMUVUVVN 


load exp 

odivisions=3 ; 

load odconf 

numniodes=8 ; 

odiength=2.001*pi ; 

firsttime=1; 

dw=odlength/(odivisions+1 ); 
dow=(freqtop-freqbottom)/(fineness-1 ); 
owref={reqbottom:dow-freqtop; 


if firsttime 
if meters 
waitbar_handle=waitbar(0,'Computing Experimental FRI"); 
else 
disp(‘Getting experimental FRF') 
end 
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H_T_EXP={]J; 
for count=1:length(owref) 
if ineters 
waitbar(count/length(owref)); 
end 
Z_exp=k_expt+)*owref(count)*c_exp-owref(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
% skyindex=1; 
% for indexl=1:aset_size 
% for index2=index|:aset_size 
% red_holder(skyindex)=h_exp(index 1 ,index2); 
% skyindex=skyindex+1; 
% end 
% end 
% H_T_EXP=[H_T_EXP red_holder']; 
H_T_EXP=[H_T_EXP h_exp(:)]; 
end 
if meters 
close(waitbar_handle) 
end 
save TPeEXPH T EXP 
clear H_T_EXP 
Yload HLANAL RED; 
Yuncorrha=H_ANAL RED; 
%clear H_ANAL RED 
uncorrha={]; 
if meters 
waitbar_handle=waitbar(0,'Computing Uncorrected FRF'); 
else 
disp('Getting Uncorrected FRF’) 
end 
load stat 
for i=1:length(owref) 
if meters 
waitbar(i/length(owref)), 
end 
tempza=kstat+j *owref(i)*cstat-owref(1)*2* mstat; 
tempha=inv(tempza); 
uncorrha=[uncorrha temphac:)]; 
end 
if meters 
close(waitbar_handle) 
end 
save uncorrha uncorrha 
clear uncorrha 


end 

load H_T_EXP 
temp=(]; 

for index=1:nunimodes 


for coord=1:aset_size 
if index == 
antiresfreqs=find(owref>0 & owref<omegax(index)); 
else 
antiresfreqs=find(owref>omegax(index-1 ) & owref<omegax(index)); 
end 
antires=min(log10(abs(H_T_EXP(ndx3d({aset_size aset_size length(owref)].... 
coord,coord,antiresfreqs))))); 
whre=find(log10(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref).... 
coord,coord,antiresfreqs))))==antires); 
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temp=owref(antiresfreqs) ; 
antiresfreq(coord,index }=temp(whre( 1 )); 
end 
end 
temp=[]; 
for index=1 :nummodes 
resfreqs=find(owref>omegax(index)-dow & owref<omegax(index}+dow), 
res=max(log10(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)].... 
cset_rel(1),cset_rel(1),resfreqs))))), 
whre=find(log10(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)].... 
cset_rel(1),cset_rel(1),resfreqs))))==res), 
temp=owref(resfreqs) ; 
resfreq(index)=temp(whre(1)); 
end 
clear H_T_EXP 


Yclear omegax 
*omegax=resfreq; 


for loopindex=startloop:skiploop:endloop 
ODZ="]; 
temp=[]; 
rdiag=[]; 
ow=[], 
lowwer=owref(1}-.5*odlength; 
uppper=owref(1)+.5*odlength ; 
% if odivisions == 
%  ow=[ow owref(1)] 
% else 
% ow=[ow lowwertdw:dw:-uppper-dw]; 
% end 
% rdiag=[rdiag ones(1,odivisions*cset_size). *owref(1)]; 


for index=1:loopindex 
%  lowwer=omegax(index).5*odlength ; 
%  uppper=omegax(index)}+.5*odlength ; 


lowwer=resfreq(index)-.5*odlength ; 
uppper=resfreq(index}+.5*odlength ; 

ow=[ow lowwertdw-.dw-uppper-dw}; 

rdiag=[rdiag ones(1 ,odivisions*cset_size).*resfreq(index)]; 
if use_antires 


for coord=1:aset_size 
lowwer=antiresfreq(coord,index)-.5 *odlength; 
uppper=antiresfreq(coord,index)}+.5*odlength; 
size(lowwer+dw.dw:uppper-dw) 
ow=[ow lowwertdw-dw:uppper-dw] 
rdiag=[rdiag ones(1,odivisions*cset_size). *resfreq(index)}; 
end 
end 
end 
R=diag(rdiag); 
load stat 
if meters 
waitbar_handle=waitbar(0,'Computing DK/DM/DC’), 
else 
disp(‘Getting DK/DM/DC’) 
end 
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200 
for count=1:length(ow) 
if meters 
waitbar(count/length(ow)); 
end 
205 z_anal_red=kstat-ow(count)*2*mstat; 


h_anal_red=inv(z_anal_red),; 
z_exp=k_exptj*ow(count)*c_exp-ow(cout)*2 *m_exp; 
h_exp=inv(Z_exp); 

h_exp=h_exp(aset,aset), 

210 hacc=h_anal_red(cset_rel,cset_rel); 
hxec=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hacc-hxcc)*inv(hacc))-hacc); 
ODZ=[ODZ; dz]; 
temp=[temp; eye(cset_size) -ow(count)*2*eye(cset_size) +j*ow(count)*eye(cset_size)]; 

215 end 

if meters 

close(waitbar_handle) 

end 

DKDMDC 1=inv(temp'*inv(R)*temp)*temp'*inv(R)*ODZ; 
220 %DKDMDC1=temp\ODZ ; 

ODK1=DKDMDC1(1:cset_size,:) ; 

ODM1=DKDMDC 1 (cset_size+1:2*cset_size,:) ; 

ODC1=DKDMDC1(2*cset_size+1:3*cset_size,:) ; 


229 clear hace hxce h_exp z_exp z_anal_red h_anal_red dz 
clear temp ODZ 


%WH%MMM%M%MM%M%N0% 09 0% 0o%o%N%2 09 0 0% 0% o%%o% 09 0% 0%%9 01%%%%% MMM Y0%? 02 0%%% 0%o%o% 09 0% V%o% o%V%WN% 
230 APWVWAHVWWYWHV%Y% 
load stat; 
odcorrha=[]; 
kstat(cset_rel,cset_rel)=kstat(cset_rel,cset_rel}}ODK1; 
mstat(cset_rel,cset_rel)}=mstat(cset_rel,cset_rel }ODM1;; 
235 cstat(cset_rel,cset_rel)=cstat(cset_rel,cset_rel}+}ODC1; 


if meters 
waitbar_handle=waitbar(0,'Computing Corrected FRF’); 
else 
240 disp(‘Getting DK/DM/DC’) 
end 
for i=1:length(owref) 
if meters 
waitbar(i/length(owref)); 
245 end 
tenpza=kstatt+j *owref(i)*cstat-owref(i)*2 *mstat; 
tempha=inv(tempza); 
odcorrha=[odcorrha tempha(:)}; 
end 
250 if meters 
close(waitbar_handle) 
end 
if loopindex == 1 
odcorrhal=odcorrha; 
255 save od] odcorrhal 
clear odcorrhal 
elseif loopindex == 2 
odcorrha2=odcorrha; 
save od2 odcorrha2 
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clear odcorrha2 


elseif loopindex == 


odcorrha3=odcorrha; 
save 0d3 odcorrha3 
clear odcorrha3 


elseif loopindex == 4 


odcorrha4=odcorrha; 
save od4 odcorrha4 
clear odcorrha4 


elseif loopindex == 


odcorrha5=odcorrha; 
save od5 odcorrhaS 
clear odcorrha5 


end 


clear tempha tempza kcorrected mcorrected ccorrected 
clear ODM1 ODC1 DKDMDC1 ODZ temp kstat mstat 


load uncorrha 
load H_T_EXP 
if complete 


plotwhichcoord=[cset(mid_index)]; 


else 


plotwhichcoord=[cset(mid_index)]; 


end 
for plotindex=1:length(plotwhichcoord) 


% 
% 
% 


which_coord=plotwhichcoord(plotindex); 
h=figure(h+1 ); 
fignum=fignum+] ; 
subplot(2,1,1); 
plot(owref/(2*p1).... 
log10(abs(H_T_EXP(nndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,"’)))),'T-,... 
owref/(2*pi).... 
log10(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,'')))),’g—,... 
owref/(2*pi),... 
log 10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)).... 
which_coord,which_coord,"’)))),’b-.") 
ylabel(('H(',int2str(aset(which_coord)).... 
' nt2str(aset(which_coord)),") in/ibf (log10 of)']) 
hold on 
v=axIS; 
Yxlabel('Omega (Hz)’) 
if titles 
if use_antires 
title({int2str(odivisions),' PT MATRIX SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 
else 
title([int2str(odivisions),’ PT MATRIX SOLTNS WEIGHTED LEFT TO RIGHT) 
end 
end 
plot(resfreq(1:loopindex)/(2*pi).... 
ones(1 ,loopindex)*(v(3 }}-abs(v(4}-v(3 ))*.98),'kt+') 
if use_antires 
%plot(diag(antiresfreq(1:loopindex, | :loopindex) (2 *pi),ones( 1 ,loopinde%x)*(v(3 }+abs(v(4)-v(3))*.98), k*') 
end 
grid on 
hh=legend(‘Experimental FRF’,'Corrected MAT Anal FRI",'Uncorrected Anal FRF',Included Modes’), 
axes(hh) 
hold off 
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Yif pswitch=="'y' 
% prtfigl(fignum) 
% delete(h) 
% M%else 
% delete(h) 
% end 
subplot(2,1,2) 
thetop=find(abs(owref-ow(length(ow))* 1.5 *ones(1,length(owref)))==min(abs(owref- 
ow(length(ow))*1.5*ones(1,length(owref))))); 


plot(owref(1 :thetop)/(2*pi).... 
log] 0(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)}.... 
which_coord,which_coord, | :thetop)))),'r-... 
owref(] :thetop)/(2*p1).... 
log 10(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,| :thetop)))),'g—"'.... 
owref( 1 :thetop)/(2*pi).... 
log10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
whuch_coord,which_coord, 1 :thetop)))),‘b-.’) 
ylabel([H(,int2str(aset(whuich_coord)),",",int2str(aset(which_coord)),") in/lbf (log10 of)'J) 
hold on 
V=anls; 
xlabel(‘(Omega (Hz)') 
if titles 
if use_antires 
title([int2str(odivisions),' PT MATRIX SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 
else 
title([int2str(odivisions),’ PT MATRIX SOLTNS WEIGHTED LEFT TO RIGHT) 
end 
end 
plot(resfreq(1 :loopindex)/(2* p1),ones(1 ,loopindex)*(v(3}+abs(v(4}-v(3 ))*.98),'k+') 
%  ifuse_antires 
% M%plot(antiresfreq(which_coord, ]:loopindex)(2*p1),ones(1 ,loopindex)*(v(%3 }+abs(v(4}-v(3))*.98),'k*") 
% end 
erid on 
%hh=legend(‘Expenmental FRF','Corrected MAT Anal FRF’,"Uncorrected Anal FRF', Included Modes’), 
Y%axes(hh) 
hold off 
if loopindex == 
print -dedjcolor 
print -dinfile fig6_1 
elseif loopindex == 
if odivisions == 
print -dcedjcolor 
if complete 
print -dmfile fig6_9 
else 
print -dnifile fig6_6 
end 
end 
elseif loopindex == 
print -dcedjcolor 
print -dinfile fig6 2 
end 
% if pswitch=='y' 
%  prtfigl(fignum) 
% delete(h) 


% else 
% delete(h) 
% end 
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380 end ' 
delete(h) 
clear H_T EXP 
clear odcorrha 
clear uncorrha 
665 end 
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WWVPCHAPE_2.MZHUVWVUVUVUVWUUVVU%VVVVUVUVYVVUVVUVUYVVVUVVVUVUwLUVVUWVHUVUVHUVUUHVU/M% 
MMW*wWVNVAHVVN% 

%%% Decompose DZ into DM, DK, & DC RwR*wR%VwVKwV%%UVVRwVUMUUNVMVUUVRwVUUUVRwVUVYUwYUwUUvVHUUUVUwUM% 
%%% using matrix formulation  %%%%V%V%VV%VVM%VwVU%KwVUVKVUVUV%VUVUUUUVVUwLUvUUVUYUHUVUWUM% 
YUWWVLYRWVUWLAPWVVLVVUVVWUVVUVUVVUVUVVLVVVLVVUVVVVUVVUVVLUVUVVwLVUUW“LUUYVUVVVUUVU% 
PPWHAVMAVVVUVUVVY 

cle 

clear 

closeall 

load INT 

whitebg(‘white') 

close 

h=0; 

fignum=1; 

use_antires=1 

titles=0 

nud_index=round(length(cset)/2); 


clear L Z_ANAL c_exp conn h_anal_red 

clear kexto m_exp mexto temp_1| diags 

clear true_damping true_mass true_stiffness z_anal z_anal_red z_exp 
clear L Z_ ANAL c_exp conn h_anal_red 

clear kexta kexto mextO temp_1_diags 


YRRVVwV’w“vVwMNMVH“V“L“MUMNMM“VM*VM“MFORCE THE CSET%H%%MMMUMMMMMMUMMMU%N% 
YRRwHV“HwMUVM“wM“N“VM*“Y%SPATIALLY INCOMPLETE%%%%%%%M%MMMMNMNNMNMNMNN% 
if complete 
YRVAwANwwWN%N’MNM“VAw“VU~N’MNNNMNMNMNMNMNMNMUMNMNMNMNMMMNMNMMMMNMMUNMMUMNMNNM%MU% 
AHYRVUwVAHwWUWNV“vV““VwHUVMV%SPATIALLY COMPLETE%N%%NMMNNNNNMNMNNMNMM%N% 
cset=[4 5 6789 10 11]%%%%%NMNMN%%NMMMMNMMMMNNMM% 
cset_rel=csethwMUMnww%%UVUMUMMAHUVMNMUVAMUUUUNMUUMVWUUMNMMV%V%% 
cset_size=length(cset)%%%%V%%%%w%UVUANMUMUUMUUUUMUMUUVAMUVUVVN% 
YRVHwVww%MNwYNMMUMnNNUM’UNUM*VMUwNMNNMNMNMNMUVNMNNMUVUNMNMVYNUNNVUUNUMUUVMNMNMUHVN 
else 
YRV“V“V“UW*wV“V*V“W%SPATIALLY INCOMPLETE%%%%%%H%NMNNMNMNMNMNMNMMM%N% 
cset=[1 3579 11 13 1ISJI*KZAwHHRHRwHUwHUVUwW%V%HNNMUHUWWUwYYYY 
cset_rel=[1 23 45 6 7 8I%NNHMMnMUMNM%WNMNNMUMUMNUNNMUNNN 
cset_size=length(cset)%%%%»%%%%%%MNMUMNNMUMUNMMNMUMKVMUUM%%% 
end 
AYRwW“wMwMNMUAwMNwMAHUUVMWwMNNMMNUMNMMNMUWMNMMNMMUMNMNMMMUMM’MNMMMUNMNMNMNMNNMN% 


PUVAWALPNLVWAVLLAWLVAAWALLLLWWVLLVVVWVLVSASWVWLVNWVL/VVVUVM% 
YWRW“W~/AWLA%%YYSPATIALLY COMPLETE%%%%%%V%%%UL%VWVWU%M% 

Yecset=[7 8 9 10 11 12 13 14]%%%%%%%%%%%W%W%~Z~YV%%~VOVO% 
Ycset_rel=cset%%VVW%VV~V%V~V~W“LULV~WWVWVWUVVWV“VVVWVU%V LVM 
Yecset_size=length(cset)~%LVV%V%WAW%V~WLVVVWLUV%V%UVVNNY% 

LUPRLLA AAV AYL ALLAN LLYVAWLVVVWUYVAVLVVV VV VV% V/V 


load exp 

odivisions=3 ; 

load odconf 

nummodes=9_ ; 

odlength=2*pi ; 

firsttime=1; 
dw=odlength/(odivisions+1 ); 
dow=(freqtop-freqbottom )/(fineness-1 ); 
owref=freqbottom:dow-freqtop; 


if firsttime 
if meters 
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waitbar_handle=waitbar(0,'Computing Experimental FRF'): 
else 

disp(‘Getting experimental FRF") 
end 


H_T_EXP={]; 
for count=1:length(owref) 
if meters 
waitbar(count/length(owref)), 
end 
z_exp=k_exptj *owref(count)*c_exp-owret(count)*2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
% skyindex=1; 
% for index1=1:aset_size 
% for index2=1ndex] :aset_size 


% red_holder(skyindex)=h_exp(index 1 ,index2), 
% skyindex=skyindex+ 1; 

% end 

% end 


%  H_T_EXP=[H_T_EXP red_holder’}; 
H_T_EXP=(H_T_EXP h_exp(:)}; 
end 


if meters 
close(waitbar_handle) 
end 


save H_T_EXPH T_EXP 
clear H_T_EXP 
end 
%load H_ANAL_RED; 
Y%uncorrha=H_ANAL_RED; 
Yclear HAANAL RED 
if firsttime 
uncorrha=[]; 
if meters 
waitbar_handle=waitbar(0,'Computing Uncorrected FRF"’); 
else 
disp('Getting Uncorrected FRF’) 
end 


load stat 
for i=1:length(owref) 
if meters 
waitbar(i/length(owref)), 
end 
tempza=kstat+j *owref(1)*cstat-owref(1)*2 *mstat; 
tempha=inv(tempza); 
uncorrha=[uncorrha tempha(:)}; 
end 


if meters 
close{waitbar_handle) 
end 


save uncorrha uncorrha 
clear uncorrha 
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eud 


load H_T_EXP 
510 for nndex=] :nununodes 
for coord=]:aset_size 
if index == 1 
antiresfreqs=find(owref>0 & owref<omegax(index)); 
else 
oS antiresfreqs=find(owref>omegax(index-1) & owref<omegax(index)): 
end 
antires=min(log]0(abs(H_T_EXP(ndx3d([1 aset_size*(aset_size+1}/2 length(owref)].... 
1 ,skyred(coord,coord),antiresfreqs))))); 
whre=find(log10(abs(H_T_EXP(ndx3d({1 aset_size*(aset_size+1)/2 length(owref)].... 
920 1 ,skyred(coord,coord),antiresfreqs))))==antires); 
temp=owref(antiresfreqs) ; 
antiresfreq(coord,index)=temp(whre), 
end 
end 
Jeo 
for index=] :numumodes 
resfreqs=find(owref>omegax(index)-dow & owref<omegax(index)}+dow); 
res=max(log]0(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)].... 
cset_rel(1]),cset_rel(1),resfreqs))))); 
530 whre=find(log10(abs(H_T_EXP(ndx3d({aset_size aset_size length (owref)].... 
cset_rel(1),cset_rel(1),resfreqs))))==res); 
temp=owref(resfreqs) ; 
resfreq(index)=temp(whre); 
end 
JOO 
clear H_T_EXP 
“clear omegax 
Yomegax=resfreq; 


540 for loopindex=startloop:skiploop:endloop 


subdivisions=9; % should be odd for simpson's rule 
intlength=2* pi; % 1 Hz bandwidth 
dW=intlength/subdivisions; % sampling freq =.25 HZ 
W=[); 

545 W1={]; 


weight=[]; 
lowwer=owref(1 }-.5*odlength; 
uppper=owref(1 }+.5*odlength; 
W=[W lowwertdW:dW:uppper-dW ]; 
550 W1=[W1 lowwert+tdW:dW:uppper-dW]; 
weight=[weight ones(size(lowwertdW:dW:uppper-dW))}; 


for index=1 :loopindex 
lowwer=omegax(index)-.5*intlength; 
999 uppper=omegax(index)+.5*intlength; 
W=[W lowwert+dW:dW:uppper-dW]; 
W1=[W1 lowwert+dW:dW:uppper-dW]; 
weight=[weight ones(size(lowwer+dW:dW:uppper-dW ))}; 
lower=antiresfreq(index)-.5*intlength; 
560 upper=antiresfreq(index)}+.5 *intlength; 
W=[W lowwert+dW:dW:uppper-dW]; 
W1=[W1 lowwert+dW:dW:uppper-dW]; 
weight=[weight oues(size(lowwertdW:dW:uppper-dW ))]; 
end 


965 
%Ycompute integrals“%%%%%M%%%MU%UVMUVV%MNMVMUMUVU%%M% MUWVUVUUVHUVUVAWVOVVYV0V0V%0V0VoV0%0 
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DZ_R=(]; 
DZ_I={]; 
if meters 
570 waitbar_handle=waitbar(0,'Computing DZ’), 
else 
disp(‘Computing DZ’), 
end 
load stat 
0/75 for count=1:length(W) 
if meters 
waitbar(count/length(W)),; 
end 
z_anal_red=kstat/(j* W(count)}+cstat/(j*W(count))+j*W(count)*mstat; 
580 h_anal_red=inv(z_anal_red); 
hacc=h_anal_red(cset_rel,cset_rel); 


Z_exp=k_exp/(j*W(count)}+c_exp/(j* W(count))}+j* W(count)*m_exp; 
h_exp=inv(zZ_exp); 

585 h_exp=h_exp(aset,aset); 
hxcc=h_exp(cset_rel,cset_rel), 


dz=inv(inv(inv(hacc)*(hace-hxec)*inv(hacc))-hacc), 
dz_r=real(dz); 
590 dz_i=imag(dz), 
DZ_R=[(DZ_R dz_1(:)}; 
DZ_I=[DZ_I dz_i(:)}; 
end 


995 if meters 
close{waitbar_handle) 

end 

W1=1./W1,; 

(INTK, INITM, INTC]=intsub(DZ_LDZ_R,W,W1,weight,cset_size), 
600 intcorrha=[]; 

clear hacc hxcc h_exp z_exp z_anal_red h_anal_red dz 

clear temp ODZ 


605 %HH}%H%%%%M%%%%%M%%%%HMMMHH%%MMMHM%MMHMHMMMM%%H%M%MM%M%MMMMMMMMMMUMMMMM%%N% 
NMW%H%%%M%MHV% 
load stat; 
odcorrha=[]; 


610 kstat(cset_rel,cset_rel)=kstat(cset_rel,cset_rel }+INTK; 
mstat(cset_rel,cset_rel)=instat(cset_rel,cset_rel +INTM, 
cstat(cset_rel,cset_rel)=cstat(cset_rel,cset_rel }+INTC; 


if meters 
615 waitbar_handle=waitbar(0,'Computing Corrected FRF’); 
else 
disp(‘Getting DK/DM/DC’) 
end 


620 for 1=1:length(owref) 
if meters 
waitbar(i/leneth(owref)), 
end 
tempza=kstat+j* owref(1)*cstat-owref(1)*2 *mstat; 
625 tempha=inv(tempza), 
odcorrha=[odcorrha tempha(:)]; 
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end 


if meters 
Close(waitbar_handle) 
end 


if loopindex == 1 
lodcorrhal =odcorrha; 
save lod| lodcorrha} 
clear lodcorrha] 

elseif loopindex == 2 
lodcorrha2=odcorrha; 
save lod2 lodcorrha2 
clear lodcorrha2 

elseif loopindex == 
lodcorrha3=odcorrha; 
save lod3 lodcorrha3 
clear lodcorrha3 

elseif loopindex == 
lodcorrha4=odcorrha; 
save lod4 lodcorrha4 
clear lodcorrha4 

elseif loopindex == 5 
lodcorrhaS=odcorrha; 
save lod5 lodcorrhaS 
clear lodcorrha5 

end 


clear tempha tempza kcorrected mcorrected ccorrected 
clear INTM INTC INTK temp kstat mstat cstat 

load uncorrha 

load H_T_EXP 


if complete 
plotwhichcoord=[cset(mud_index)]; 
else 
plotwhichcoord=[cset(nud_index)}; 
end 


for plotindex=1 :length(plotwhichcoord) 

which_coord=plotwhichcoord(plotindex); 

h=figure(h+1 ); 

subplot(2,1,1) 

plot(owref/(2*p1).... 
logl0(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,""’)))),'r-‘,... 

owref/(2*p1),... 
log10(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,')))),'g—... 

owref/(2*pi),... 
log10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,":’)))),’b-.') 

ylabel((H(,int2str(aset(which_coord)),... 
’, mnt2str(aset(which_coord)),") in/lbf (log10 of)'}) 

hold on 


V=axis; 
Yxlabel(‘Omega (Hz)') 
if titles 


if use_antires =='on' 
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an (int2str(odivisions),' PT INTEGRAL SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 
else 
title([int2str(odivisions),’ PT INTEGRAL SOLTNS WEIGHTED LEFT TO RIGHT)) 

end 
end 
plot(omegax(1:loopindex)/(2*pi).... 

ones(1,loopindex)*(v(3}+abs(v(4}-v(3))*.98),'k+') 

if use_antires 

“plot(antiresfreq( which_coord,]:loopindex)/(2*pi),ones(1,loopindex)*(v(%3 )+abs(v(4)-v(3))*.98),'k*') 

end 
grid on 
hh=legend(‘Experimental FRF",'Corrected INT Anal FRF',Uncorrected Anal FRF’,lucluded Modes’), 
axes(hh) 
hold off 
Yif pswitch=="'y' 
% prtfigl(fignum) 
% delete(h) 
%else 
% delete(h) 
%end 
subplot(2,1,2) 
thetop=find(abs(owref-W(length(W))* 1 .3*ones(1 ,length(owref)) )}==min(abs(owref- 


Wlength(W))* 1.3*ones(1,length(owref))))); 


% 
% 
% 


% h=figure(h+] ); 
% fignum=fignumt+| ; 
plot(owref(1 :thetop)/(2*pi).... 
log] 0(abs(H_T_EXP(ndx3d({aset_size aset_size length(owret)].... 
which_coord,which_coord,1:thetop)))),'r-’.... 
owref(1:thetop)/(2*pt).... 
log! 0(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
Wluch_coord,which_coord,1:thetop)))),’g—"’.... 
owref(1:thetop)/(2*pi).... 
log10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord, 1:thetop)))),’b-.’) 
ylabel({H(‘int2str(aset(which_coord)),',’int2str(aset(which_coord)),') uv/lbf (log10 of)']) 
hold on 
V=axis; 
xlabel((Omega (Hz)’) 
if titles 
if use_antires == ‘on ' 
title([int2str(odivisions),’ PT INTEGRAL SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 
else 
title((int2str(odivisions),’ PT INTEGRAL SOLTNS WEIGHTED LEFT TO RIGHT) 
end 
end 
plot(omegax(1 :loopindex)/(2* p1),ones(1,loopindex)*(v(3)+abs(v(4}-v(3))*.98),'k+’) 
if use_antires 
%plot(antiresfreq( which_coord, ]:loopindex)/(2* pi),ones(1,loopindex)*(v(%3 }-abs( v(4 )-v(3))*.98),'k*") 
end 
grid on 
%hh=legend(‘Experimental FRF','Corrected INT Anal %RF',Uncorrected Anal FRF',Iucluded Modes’); 
%Yaxes(hh) 


if loopindex == 
print -dcdjcolor 
print -dinfile fig6_3 
h=figure(h+1 ); 
hold off 
load od2 
plot(owref(1:thetopV(2* pt)... 
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logl0(abs(H_T_EXP(ndx3d({aset_size aset_size length(owref)].... 
which_coord,which_coord,1:thetop)))),'r-.... 
owref(1:thetop)/(2*p1).... 
log10(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord, | :thetop)))),'g—,... 
owref( 1:thetop)/(2*pi).... 
log] 0(abs(odcorrha2(ndx3d([aset_size aset_size length(owref)).... 
which_coord,which_coord,1:thetop)))),’b-."... 
owref(1:thetop)/(2*pi).... 
log10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,1:thetop)))),'m:') 
ylabel(('H(‘,int2str(aset(which_coord)),',’,int2str(aset(which_coord)),") in/Ibf (log10 of)'}) 
hold on 
V=axis, 
xlabel('Omega (Hz)') 
if titles 
if use_antires == ‘on ' 
title({int2str(odivisions),' PT INTEGRAL SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 
else 
title({int2str(odivisions),’ PT INTEGRAL SOLTNS WEIGHTED LEFT TO RIGHT) 
end 
end 
plot(omegax(1:loopuidex)/(2*pi),ones(1,loopindex)*(v(3)}+abs(v(4}-v(3))*.98), K+) 
% if use_antires 
% Yplot(antiresfreq(which_coord, 1 :loopindex)/(2*pi),ones(1 ,loopindex)*(v(%3 }tabs(v(4}-v(3 ))*.98), k*') 
% end 
grid on 
hh=legend(‘Experimental FRF','Corrected INT Anal FRF','Corrected MAT Anal FRF',Uncorrected Anal FRF', ‘Included 
Modes’); 
axes(hh) 
hold off 
print -dedjcolor 
print -dmfile fig6_5 
elseif loopindex == 3 
if odivisions == 1 
print -dedjcolor 
print -dmnifile fig6_7 
h=figure(h+1); 
hold off 
load od3 
plot(owref(1:thetop)/(2*pi).... 
log10(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)).... 
which_coord,which_coord, 1:thetop)))),'t-’,... 
owref(1:thetop)/(2*p1).... 
log10(abs(odcorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,]:thetop)))),'g—\... 
owref(1:thetop)/(2*p1).... 
log10(abs(odcorrha3(ndx3d([aset_size aset_size length(owref)].... 
which _coord,which_coord,1:thetop)))),b-.’,... 
owref(1:thetop)/(2* pi)... 
log 10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
Which_coord,which_coord, 1:thetop)))),'m:') 
ylabel(('H(‘,int2str(aset(which_coord)),’,',int2str(aset(which_coord)),') uvIbf (log!0 of)'}) 
hold on 
v=anis, 
xlabel(‘Omega (Hz)') 
if titles 
if use_antires == ‘on' 
title({int2str(odivisions),' PT INTEGRAL SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT) 


else 
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title([int2str(odivisions),' PT INTEGRAL SOLTNS WEIGHTED LEFT TO RIGHT) 
end 


end 
plot(omegax(1 :loopindex)/(2 *pi),ones( 1 ,loopindex)*(v(3 }tabs(v(4)-v(3))*.98),'k+') 
% if use_antires 
%N “plot(antiresfreq(which_coord,1:loopindex)/(2*pi),ones(1 ,loopindex)*(v(%3)+abs(v(4)-v(3))*.98),'k*") 
% end 
grid on 
hh=legend(‘Expermmental FRF','Corrected INT Anal FRF',"Uncorrected Anal FRI"""Included Modes’), 
axes(hh) 
hold off 
print -dcdjcolor 
pnnt -dmfile fig6_ 8 


end 
elseif loopindex == 
print -dcdjcolor 
print -dmfile fig6_4 
end 
% if pswitch=='y' 
% prtfig 1(fignum) 
% delete(h) 
% else 
% delete(h) 
% end 
end 
delete(h) 
clear H_T_EXP 
clear odcorrha 
clear uncorrha 
end 


load H_T_EXP 
load od] 
load od2 
load od3 
load od4 
load uncorrha 
%load lod5 
h=figure(h+1); 
plot(owref/(2*p1).... 
log10(abs(H_T_EXP(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,':')))),'r-',... 
owref/(2*p1).... 
log 10(abs(odcorrhal (ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,"')))),'g—',... 
owref/(2* pi)... 
log 10(abs(odcorrha2(ndx3d([aset_size aset_size length(owref)}.... 
which_coord,which_coord,':')))),'b-..... 
owref/(2*pi).... 
log 1 0(abs(odcorrha3(ndx3d([aset_size aset_size length(owref)].... 
Wwhich_coord,which_coord,":')))),'m—".... 
owref/(2* pt)... 
log10(abs(odcorrha4(ndx3d([aset_size aset_size length(owref)],... 
which_coord,which_coord,":')))),'c-.',... 
owtef/(2*pi).... 
log 10(abs(uncorrha(ndx3d([aset_size aset_size length(owref)].... 
which_coord,which_coord,":')))),k-.') 
ylabel({HC,int2str(aset(which_coord)).... 
', int2str(aset(which_coord)),") un/lbf (log10 of)']) 
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hold on 
v=anis; 
Yxlabel(‘Omega (Hz)') 
if titles 
if use_antires == 'on' 
title({int2str(odivisions),’ PT INTEGRAL SOLTNS WITH ANTIRESONANCES WEIGHTED LEFT TO RIGHT)) 
else 
title{[int2str(odivisions),’ PT INTEGRAL SOLTNS WEIGHTED LEFT TO RIGHT) 
end 
end 
% plot(omegax(1:loopindex)/(2*p1).... 
%  ones(1,loopindex)*(v(3)}+abs(v(4)-v(3))*.98),'k+') 
% if use_antires 
% Yplot(antiresfreq(which_coord, | :loopindex)/(2*p1),ones(1 ,loopindex)*(v(%3 }}-abs(v(4)}-v(3))*.98),'k*") 
end 
V=ax1S; 
axis({10 300 v(3) v(4)]); 
grid on 
hh=legend(‘Experimental FRF","1 mode MAT solutiou','2 modé MAT solution’,'3 inode MAT solution’,'4 mode MAT 


solution’, Uncorrected Anal FRF'); 


axes(hh) 

hold off 

%if pswitch== 

% prtfigl(fignum) 
% delete(h) 
“*else 

% delete(h) 
%end 


print -dedjcolor 
print -dinfile fig6_12 
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YUASA AVVV LLLP OWA SIA AVA VAS V VASLUI VoVVVOYO%O%%% 
LYVA YMA AAA 


% PROGRAM: SETUP.M % 
% INITILIZES DATA FILE FOR A FINITE ELEMENT ANALYSIS % 
% VARIABLES SAVED TO FILE SETUP.MAT % 


PP APYPPPYVVPVHLAHWHWLVYLLLUWAHVRHUVVVVWV“LUYVUWAVVVVVHUVYUUVHVVVUUVLVUVUVVVUVU%M% 
P~YPRYYAHRVY% 


cle 
clear 
L=input(Enter the total length of the beam: '); 
numel=input(‘Num of elements: '); 
dof node=input(‘num dof per node: '), 
area=input(‘area of beam: '), 
eell=input(EI of beam: '); 
%ee=input(‘modulus for beam: '); 
pho=input(‘mass density for beam: '); 
numdof=dof_node*(numel+1) 
cle 
whule 1 
conforce=input(Enter # of concentrated loads: "); 
if ~isnan(conforce) 
break 
end 
end 
for 1=1:conforce 
forcepos(i)=input([‘Beam position of force ',num2str(1),' : ‘}); 
forcesiz(i)=imput((Magnitude of force ‘;num2str(1)," : ']); 
end 
force=zeros(numdof, 1 ); 
for 1=1:conforce 
temppos=ceil(forcepos(i)/(L/numel)); 
force(dof_node*temppos+1 )=force(dof_node*temppos+ 1 }+.5*forcesiz(1); 
force(dof_node*temppost+3)=force(dof_node*temppos+3}+.5*forcesiz(i); 


end 
while 1 
lmass=input(‘Enter # of Lumped masses: '); 
if ~isnan(lmass) 
break 
end 
end 
for 1=1:lmass 
masspos(i)=input((‘Beam position of Mass ',num2str(i),’ : '}); 
masssiz(i)}=input([‘Magnitude of Mass 'snum2str(1)," : ']); 
end 
lumpmass=zeros(numel, | ); 
for 1=1:lmass 
temppos=round(masspos(1)/(L/numel)), 
lumpmass(temppos)=lumpmass(temppos)+masssiz(1); 
end 
lumpmass' 
while 1 
Ispring=input(Enter # of Lumped Springs: ‘), 
if ~isnan(Ispring) 
break 
end 
end 
for i=1:lspring 
springpos(1)=input([‘Beam position of Spring ';num2str(1),’ : ']); 
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springsiz(1)=iput(['Spring constant of Spring ';num2str(1),' : ']); 
end 
lumpspring=zeros(numiel, | ); 
for i=1:lspring 
temppos=round(springpos(1)/(L/numel)); 
lumpspring(temppos)=lumpspring(temppos )}+springsiz(1); 
end 
lumpspning' 
conn=[1,2]; 
for 1=2:numel 
conn=[conn;i,i+1}; 


end 
com 
while 1 
bc=['pinned-pinned_' 
‘clamp-clamp—s' 
‘left guided clamp ' 
‘right guided clamp’ 
‘cantilevered ' 
‘free-free ds 
cle 
help betext 
n=input(‘Select a boundary condition: ‘); 
if ((n > 0) & (n < 7)) 
break 
end 
end 
ifn == 
be='pp’; 
elseif n == 
='cc’; 
elseif n == 3 
be='Ic’; 
elseif n == 
be='re'}; 
elseif n== 
bc="'cl' 
else 
be="fF; 
end 
clear 1 
clear temppos 
clear conforce 


clear forcepos 
clear forcesiz 
Clear Imass 
clear masspos 
clear masssiz 
clear Ispring 
clear springpos 
clear springsiz 
save setup.mat 
x=0:L/numel:L; 
for i=] :length(x) 
hai)=5; 
end 
cle 
hold off 
%Yaxis(‘off) 
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“oplot((],[]) 
Y%hold on 


plot(x,h,x,h,'x’') 
125 
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MAYYVWVVV%V%V%%%%%%VMVMV%VUVUVVVMVUVV%M%%%MMNM%MVMUAUMUUVV%VUUUVwHVUYUVVUHUYUVUUV%V% 
UWV%V%UW%V%U% 
*FORM FE BEAM MODEL WITH ERRORS AT PRESCRIBED LOCATIONS%%%%%%%%%%% 
%MASS ERROR LOCATION SPECIFIED VIA POSOFMASSER. VALUE OF ERROR IN%% 
%PERCENT GIVEN BY VALOFMASSERR, ETC.%%%%%%%%MUUUUVMUVUNMUUMUVUVNUVVUMVUV%%% 
WNWwWwW%W%%V/HVHV%W%M%M%%VWVMM%MAMM%UW%M%V%%M%M%%M%MUMUMUVVU%NM%UMUMUMMNMNMUVMUMUMUM%MVMUVUMUVMU%M% 
RWWVUVVMVY% 
clear,cle;clg 
casename="Spatially Incomplete A set Stiffness error’; 
if exist(‘setup.mat’) 
load setup Y%load existing input data file 
else 
setup “Acreate new input data file 
end 
struc_damping=0.00000000001; 
lumpdamp=lumpspring*0 : 
posofMasserr=[3 4] ; 
valueofMasserr=[0.25 0.25]; 
posotStifferr=[3 4]; 
valueofStifferr=[0.25 0.25] ; 
posofDamperr=[3 4]; 
valueof[Damperr=[0.25 0.25]; 
for i=] :length(posofMasserr) 
lumpmass(posofMasserr(i))=valueofMasserr(1), 
end 
for i=1:length(posofStifferr) 
lumpspring(posofStifferr(i))=valueofStifferr(i); 
end 
for i=1:length(posofDamperr) 
lumpdamp(posofDamperr(i))=valueofDamperr(1); 
end 
elen=L/numel; 
if bc=="'pp' 
doftokill=2; 
elseif be=='cc' 
doftokill=4; 
elseif bc==' 
doftokill=0; 
elseif be=='cl' 
doftokill=2; 
end 


YVAN AYA AYALA LAVAL AVAL 
%BUILD THE ELEMENTAL MASS AND STIFFNESS MATRICES%%%%%%% 

%FOR A 2DOF/NODE "F.E." STRUCTURES. %%%%%%%%~%V%W%V/%*LL%%%V0% 

WLUW AVVYV AVY ALLAN WLLL VALLLLVVVVAV VV WV%V%Y 


“%ELEMENT STIFFNESS MATRIX%%%%%%H2%%%MVMNMUVUMVUANMUMUMUMVNMVYUANNM%N% 
ke=[12 6*elen -12 6*elen; 
6*elen 4*(elen’2) -6*elen 2*(elen’2),; 
-12 -6*elen 12 -6*elen; 
6*elen 2*(elen*2) -6*elen 4*(elen’2)]}; 


ke=(eeii/(elen’3)).*ke; 
g=386; 
elemke=eeii/(elen’3); 


%ELEMENT MASS MATRICE%%%%%%%%M%M%NM%V%MMVYVMVVYVYYYOY% 


me=[156 22*elen 54 =l3-elen: 
22*elen 4*(elen’2) 13*elen -3*elen“%2; 
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54 + 13*elen 156 = -22* len; 
-13*elen -3*(elen’2) -22*elen 4*(elen’2)]; 


ime = (((pho*area)* elen/g)/(420)).*me; 
elemme=(((pho* area)* elen/g)/(420)), 


YRWYRwRHwW’>’wVHMVwMUHVUMU*wWUMUWUVMU*AHUWUVVUV%MNMMUVMVMNMNNMAMUMUMVMAHVWUMMUWMUVMUVVNUNUAVZAV%V%AVWMUAVUVV%% 
SSA SANA SASK SASASASASA: 
%ASSEMBLE GOBAL STIFFNESS AND MASS MATRIX FOR A 2 DOF F.E. STRUCTURE%%%% 
%BASED ON THE ELEMENTAL MATRIXES KE AND ME. THE STRUCTURE CONSISTS OF%%% 
%NUMEL ELEMENTS WITH ELEMENT CONNECTIVITY GIVEN BY MATRIX CONN%%%%%%%%%% 
Y~WVUVYVwHVHVYVV%VMUWM%%V%HVUVVWVMUVVUWMUUVV%MUVVMUVMUUUMVWUHV~VVAMMVMVUYVYUVUVVUVUVYUVYVUUVVAV%% 
NLA SASL SASL NASA SA SAYA: 
goblk=zeros(numdof); 
goblin=zeros(nuimdof), 
goblce=zeros(numdof); 
for i=] :numel 
v=conni(1, | ); 
w=conn(i,2); 
goblk(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v}-... 
goblk(dof_node*v-1:dof_node*v,dof_node* v-1:dof_node*v)-... 
ke(1:dof_node, | :dof_node); 
goblk(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)-... 
goblk(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)-... 
ke(1:dof_node,dof_nodet+1:2*dof_node); 


goblk(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)-... 
goblk(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)-... 
ke(dof_nodet1:2*dof_node, l:dof_node); 
goblk(dof_node*w-1:dof_node* w,dof_node*w-]:dof_node*w)-... 
goblk(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w).... 
ke(dof_node+1:2*dof_node,dof_nodet+1:2*dof_node); 


goblin(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v)}=... 
goblin(dof_node*v-1:dof_node*v,dof_node* v-1:dof_node*v)-... 
me(1:dof_node, 1 :dof_node); 
goblm(dof_node* v-1:dof_node*v,dof_node*w-|:dof_node*w)-... 
goblm(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)... 
me(1:dof_node,dof_node+] :2*dof_node),; 
goblm(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)=... 
goblm(dof_node*w-|:dof_node*w,dof_node* v-1:dof_node*v}.... 
me(dof_nodet1:2*dof_node, 1:dof_node), 
goblin(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w}... 
goblm(dof_node*w-1:dof_node*w,dof_node* w-1:dot_node*w)-... 
me(dof_nodet+1:2*dof_node,dof_node+1:2*dof_node), 

end 

goble=sqrt(-1)*struc_damping.* goblk; 

goblkx=goblk,; 

goblcx=goblc; 

goblinx=goblin; 


for i=] :numel 
v=conn(i, 1); 
W=coi(1,2); 
goblkx(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v)=... 
eoblkx(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v}-... 
jumpspring(i). *ke(1:dof_node,1:dof_node), 


goblkx(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)=... 
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goblkx(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)-... 
lumpspring(1).*ke(1:dof_node,dof_node+1:2*dof_node); 


goblkx(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)=... 
goblkx(dof_node*w-|:dof_node*w,dof_node*v-1:dof_node*v)... 
lumpspring(1).*ke(dof_node+1:2*dof_node,1:dof_node); 


goblkx(dof_node*w-1:dof_node*w,dof_node*w-|I:dof_node*w)-... 
goblkx(dof_node*w-1|:dof_node*w,dof_node*w-1:dof_node*w)-... 
lumpspring(1).*ke(dof_nodet+1:2*dof_node,dof_node+1:2*dof_node), 


goblmx(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v)=... 
goblmx(dof_node*v-1:dof_node*v,dof_node*v-1]:dof_node*v}... 
lumpmass(1). *me(1:dof_node, l:dof_node); 


goblmx(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)-... 
goblmx(dof_node*v-1]:dof_node*v,dof_node*w-1:dof_node*w)-... 
lumpmass(1).*me(1:dof_node,dof_node+1:2*dof_node); 


goblmx(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)=... 
goblmx(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)t... 
lumpmass(i).*me(dof_nodet+1:2*dof_node,1:dof_node), 


goblmx(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w)-... 
goblmx(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w)-... 
lumpmass(1).*me(dof_node+1:2*dof_node,dof_node+1:2*dof_node); 


goblcx(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v)=... 
goblcx(dof_node*v-1:dof_node*v,dof_node*v-1:dof_node*v)+... 
lumpdamp(i)*sqrt(-1)*struc_damping. *ke(1:dof_node, l:dof_node); 


goblcx(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)=... 
goblcx(dof_node*v-1:dof_node*v,dof_node*w-1:dof_node*w)-... 
lumpdamp(i)*sqrt(-1)*struc_damping. *ke(1:dof_node,dof_node+1:2*dof_node); 


goblcx(dof_node*w-|:dof_node*w,dof_node*v-1:dof_node*v)=... 

goblcx(dof_node*w-1:dof_node*w,dof_node*v-1:dof_node*v)t... 

lumpdamp(i)*sqrt(-1)*struc_damping.*ke(dof_nodet1:2*dof_node,1:dof_node); 

goblcx(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w)=... 

goblcx(dof_node*w-1:dof_node*w,dof_node*w-1:dof_node*w)t... 

lumpdamp(i)*sqrt(-1)*struc_damping.*ke(dof_node+1:2*dof_node,dof_node+1:2*dof_node); 
end 


numdof=nwmdof-doftokill,; 
[gk,gm,gc,k_anal,m_anal,c_anal]=fixbces(goblk,goblm, goblc, bc), 


[ekx,emx,gcx,k_exp,in_exp,c_exp]}=fixbes(goblkx,goblmx,goblcx,bc); 
save beamdata 
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function [kstat,mstat]=fstatic(k,m,oset,aset) 
aset_size=length(aset); 

kaa=k(aset,aset); 

kao=k(aset,oset); 

koo=k(oset,oset); 

koa=kao’, 

clear k; 

k={kaa,kao;koa,koo]; 


maa=m(aset,aset), 
1ao=m(aset,oset); 
moo=m(oset,oset); 
moa=mao" 

clear m; 
=[maa,mao;moa,moo]; 


t_static=-koo\koa; 
T_static=[eye(aset_size); t_static]; 


kstat=T_static'*k*T_static; 


_ mstat=T_static'*m*T_static; 


end 
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% 

function (kirs,mirs]=firs_tam(k,m,oset,aset) 

% 

% this function returns the IRS reduced stiffness 
% and mass matrices, given the unreduced couterparts. 
% Care must be taken that the aset and oset vectors correspond 
% with the existing arrangement of k and m. 

% k and m are UNPARTITIONED matrices. 

% 

aset_size=length(aset); 

% 

kaa=k(aset,aset), 

kao=k(aset,oset); 

koo=k(oset,oset); 

koa=kao’; 

clear k; 

k=[koo,koa;kao,kaa]; 

% 

maa=ni(aset,aset); 

niao=mi(aset,oset); 

moo=m(oset,oset); 

moa=mao*, 

clear m, 

m=[moo,mioa;niao,maa]; 

% 

t_static=-koo\koa; 

T_static = [t_static; eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 
mstat=T_static'*m*T_static; 

% 
tirs=t_statictinv(koo)*(moa+moo*t_static)*iv(mstat)*kstat; 
T_irs=(tirs;eye(aset_size)]; 

% , 

kirs=T_irs'*k*T_irs; 

murs=T_irs'*m*T_irs; 

% 

% end function firs_tam 
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YWW“MUwMVMVUWMUVUHUWV“LKVMUWMUWMUHUVUMNUUMNMM%MMMUWNM%VMUVMNMMNMMNUMUWVUMNMUVwWUVUWUUMRwWUVUVUVHUMUVUMUUVNUUV%% 
YPUWRWKW“YVYKHUVVWVN% 

%this function returns a vector U containing modal%%w%w%%%%%%%%%%wNM%UNMUVUUNM%% 

%frequencies (rad/sec)*2 in ascending order along%%%%%%%%%%%%%%UZ%%%%%%%% 

%with associated mode shapes%~%%%%%V%UUM%NMUMKVUHUMUM%MNMUMNMUMNUMNHUVAVM“LUVVUUVMUUVUMNVMUNM“VU% 
AWWwKhH“H“MNWNVUWHVU“VUHVNVUWNWUVHwVUVUHNYVMUNVUHLHVAHUWVUWUVAHHVA“VVMUHVVMVYUV“V“VU~UNUHUMUVNMUVUVUVUVVN% 
RW~W“Z“MNWMKVUVUV%%N% 


function [ul ,lambda,index]=freqmode(k,m) 
(u,lambdat]=eig(m\k); 
(lambda, index]=sort(diag(lambdat)); 
ll=zeros(length(k),length(k)); 
for 1=1:length(k); 
11(1,1)=lambda(1); 
end 
lambda=diag(Il1); 
for 3 = 1:length(k) 
ul(:.j) = u(,index(j)); 
end 
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NDX3D.M 


function [r,c] = ndx3d(siz,i,j,k) 


%NDX3D Index into 3-D matrix packed in a 2-D matrix. 

% ELEM = NDX3D({M N P],LJ,K) retums the element position ELEM 
% of the (i,j,k) elements of a M-by-N-by-P matrix which is 

% stored in a (M*N)-by-P matrix. For example, the three M-by-N 
% matrices Al,A2,A3 are packed into a 2-D matnx using 

% A = [Al1(:) A2(:) A3(:)]; 

% If length(I) is m, LENGTH(J) is n, and LENGTH(K) is p, 
% then ELEM will be an im-by-n-by-p matrix. 

% 

% [R,C}] = NDX3D((M N P],1,J,K) retums the row and column 
% position of the (i,j,k) element as stored in the normal 

% 2-D matrix of size (M*N)-by-P. 

% 

% To specify all the elements along one dimension use “’. 

% For instance, NDX3D([3 5 4],2:3,'7,3:4) returns the 

% elements for the 2-by-5-by-2 matrix. 

% 

% See also ELEM3D, MESHGRID, SLICE. 

% Clay M. Thompson 1 1-3-92 

% Copynght (c) 1992 by The MathWorks, Inc. 

% $Revision: 1.7$ $Date: 1993/09/03 14:36:52 § 


if isstr(1), 1 = 1:siz(1); end 
if isstr(j), j = 1:s1z(2); end 
if isstr(k), k = 1:s1z(3); end 


if isempty(i) | isenrpty(j) | isempty(k), r = []; ¢ = []; return, end 


if any( (i<O) | (i>siz(1)) ), error(‘Index I out of range.'), end 
if any( (j<0) | (j>s1z(2)) ), error(‘Index J out of range.’); end 
if any( (k<0) | (k>siz(3)) ), error(‘Index K out of range.'); end 


if nargout==2 
{jj.1i] = meshgrid(j,1), 
r = 11(:) + Gjj(:-1)*s12z(1); 
c= k(:); 
else 
{jj .11,kk] = meshgrid(j.i,k); 
r= ii + (jj-] )*siz(1) + (kk-1)*prod(siz(1:2)), 


end 
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INTSUB.M 


fiuiction [K,M,C]=intsub(Z_I, Z_R, omega, omegal, W, set_size) 
waitbar_handle=waitbar(0,'Computing Integrals’), 
jJ=sqrt(-1); 


integ=omega.*2.*abs(W), 
rint=real(integ), 
lint=1mag(integ); 
aa=mytrapz(omegal ,rint}+mytrapz(omegal ,iint)*j; 
integ=(1./omega.*2).*abs(W),; 
nut=real(integ ); 
lint=1mag(integ); 
bb=mytrapz(omegal ,rint}+mytrapz(omegal ,int)*)j; 
integ=abs(W); 
rint=real(integ); 
lint=imag(integ ); 
cc=mytrapz(omegal ,rint)}+mytrapz(omegal ,1nt)*); 
for 1=1:set_size 
waitbar(1/set_size); 
for k=i:set_size 
Mc,k)=(1/(aa* bb-ce%2))*(bb* mytrapz(omegal,... 
omega.*Z_I(ndx3d([set_size set_size length(omega)].... 
1,k, 1 :length(omega))).*abs(W))-cc*mytrapz(omegal.... 
1./omega.*Z_I(ndx3d([set_size set_size length(omega)... 
1,k, 1 :length(omega))).*abs(W))) ; 
M(k,i)=M(i,k); 
C(i,k}-(1/cc)*mytrapz(omegal,Z_R¢... 
ndx3d([set_size set_size length(omega).... 
1,k, 1 :length(omega))). *abs(W)); 
C(k,1)=C(i,k); 
Ka,k)=(1/(aa*bb-ce*2))*(cc*mytrapz(omega.,... 
omega.*Z_I(ndx3d([Sset_size set_size length(omega)]... 
1,k, 1 :length(omega))).*abs(W))-aa*mytrapz(oniega ].... 
]./omega.*Z_I(ndx3d([set_size set_size length(omega)]... 
1k, 1 :length(omega))).*abs(W))) ; 
K(k,1)=K(.k); 
end 
end 
close( waitbar_handle) 
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MYTRAPZ.M 


function z = trapz(x,y) 

%TRAPZ Trapezoidal numerical integration. 

%  Z=TRAPZ(X,Y) computes the integral of Y with respect to X using 
% trapezoidal integration. XK and Y must be vectors of the same length, 
% or X must be acolumn vector and Y a matrix with as many rows as X. 
% TRAPZ computes the integral of each column of Y separately. 

% The resulting Z is a scalar or a row vector. 

% 

%  Z=TRAPZ(Y) computes the trapezoidal integral of Y assuming unit 
% spacing between the data points. To comipute the integral for 

% spacing different from one, multiply Z by the spacing increment. 

% 

% See also SUM, CUMSUM. 


% Clay M. Thompson, 10/16/90; Cleve Moler, 1/19/92. 
% Copyright (c) 1984-94 by The MathWorks, Inc. 


% Make sure x and y are column vectors, or y is a matrix. 


% Trapezoid sum computed with vector-matrix multiply. 
{m,n]J=size(x) ; 
z=0 x 
for index=1:m 
yy=y(1+(index-1)*n:n+(index-1)*n); 
ime) ): 
Xx=x(index,:); 
XX = xx(:); 
z =z+diff(xx)' *(yy(1:n-1) + yy(2:n)y/2; 
end 
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